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ABSTRACT 



Aims. We investigate how unified models should be built to be able to predict the matter-density bispectrum (and 
power spectrum) from very large to small scales and that are at the same time consistent with perturbation theory at 
low k and with halo models at high k. 

Methods. We use a Lagrangian framework to decompose the bispectrum into "3-halo", "2-halo", and "1-halo" contribu- 
tions, related to "perturbative" and "non-perturbative" terms. We describe a simple implementation of this approach 
and present a detailed comparison with numerical simulations. 

Results. We show that the 1-halo and 2-halo contributions contain counterterms that ensure their decay at low k, as 
required by physical constraints, and allow a better match to simulations. Contrary to the power spectrum, the standard 
1-loop perturbation theory can be used for the perturbative 3-halo contribution because it does not grow too fast at 
high k. Moreover, it is much simpler and more accurate than two resummation schemes investigated in this paper. We 
obtain a good agreement with numerical simulations on both large and small scales, but the transition scales are poorly 
described by the simplest implementation. This cannot be amended by simple modifications to the halo parameters, 
but we show how it can be corrected for the power spectrum and the bispectrum through a simple interpolation scheme 
that is restricted to this intermediate regime. Then, we reach an accuracy on the order of 10% on mildly and highly 
nonlinear scales, while an accuracy on the order of 1% is obtained on larger weakly nonlinear scales. This also holds 
for the real-space two-point correlation function. 

Key words, gravitation; cosmology: theory - large-scale structure of Universe 



1. Introduction 

According to standard cosmological scenarios, the large- 
scale structures of the present Universe have formed 
through the amplification by gravitational instabil- 
ity of small almost-Gaussian primordial fluctuations 
(| Peebles! |1980[ ). Then, from obser vations of the recent 
Universe, throu gh galaxy surveys (iTegmark et al. 2006 : 



Verde et all l2002t iTakada fc JainI 120031: iKavo et al 



Scoccimarro &: CouchmanI 20011: Bernardeau et al. 2002bl: 



2004 'Sefusat ti etal 
2007t [Nishimichi et al.l 
Nishimichi et al.ll2010[) . 



2007[ 
12007 : 



Sefusatti fc Komatsu 
ISefusatti et al.l 120101: 



Cole et al."'2005'), weak-lensing studies (jMassev et al.]l2007 
Munsh i ct a l. 2008), measures of baryon acoustic oscilla- 
tions (jEisenstein et al.l 119981 |2005[) , one can derive con- 
straints on the cosmological parameters (such as the mean 
matter and dark energy contents) and on the properties 
of the initial perturbations. The main statistical quantity 
used in this context is the density two-point correlation 
function, or power spectrum in Fourier space. This allows 
measuring the scale dependence and the amplitude of the 
initial conditions, and for Gaussian fields this fully deter- 
mines all statistical properties of the matter distribution. 
However, to constrain the possible non-Gaussianities of the 
initial perturbations, to check the gravitational clustering 
scenario (through the measure of the non-Gaussianities 
generated by the nonlinear dynamics), and to break 
parameter degeneracies, it is useful to study higher order 
statistics. The three-point correlation function or bispec- 
trum in Fourier space, which is the lowest order statistics 
beyond t he Gaussian, i s the main quantity studied i n this 
context (fPeeblcs. 1980l: iFrieman fc Gaztanaga 1994 Il999t 



On large scales, or at early times, where the den- 
sity fluctuations a re smal l with in cold dark matter 
(CDM) scenarios (iPeebled [1981 . one can use per- 
turbation theory teoroff et al.l 19861: Scoccimarro! 119971 : 
^Scoccimarro ct al. 1998; Ber nardeau et al.! l2002a|) In or- 
der to extend the validity of perturbative predic- 
tions to somewhat smaller scales and to increase their 
theoretical accuracy, ma ny resummation schemes ha ve 
been proposed r ecently (j Croccc fc Scoccimarro' '2006b''a ; 
Valageas 2007a: Matarrese fc P ictroni 2007; Pictroni 200l 
Taruya fc Hiramatsull2008l : lMatsubarall2008b!:lTaruva et^ 
i2009ll . Although these methods follow different routes and 
involve different expansion schemes, they are consistent 
with each other and with the standard perturbation the- 
ory up to the order of truncation, and only differ by higher 
order contributions. Thus, they complete the standard per- 
turbation theory contributions (which consist of a finite 
number of Feynman diagrams) by (usually infinite) par- 
tial resummations of higher order diagrams. Most of these 
stu dies have f ocused on the matter power spectrum, except 
for IValageaj (|2008D , who a lso considers three-point (and 
higher) density correlations. [Bernardeau et al.l ([20081 ). who 
study higher order response functions (propagators) and 
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the bispectrum, and iBartolo et aLl (|2010[ ). who consider 
both the power spectrum and bispectrum for Gaussian and 
non-Gaussian initial conditions. 



On small scales, where the density fluctuations are 
large, one must use numerical simulations or phenomeno- 
logical models, such a s th e Lagrangian mapping in- 
troduced in Hamilton et al.l (11991 ) or the halo model 
(iMcClell and & Silk 19771: IScherrer fc Bertschinged 119911: 
iCoorav fc Sheth.2002ir Both approaches have already been 
used for the bispectrum. However, while the first method 
has only h ad som e success on intermediate mildly nonlin- 
ear scales ()Pan et a h. 2007, ), the second rr iethod is also able 
to describe the highly nonlinear scales (iMa fc Frvir2000bl: 
Scoccimarro et all boOlt IWang et al] 120041: iFosalba et al.l 



2005E 



In order to compare theoretical predictions with obser- 
vations, which go from large linear scales to small highly 
nonlinear scales, it is useful to build unified models, which 
combine for instance perturbation theories with halo mod- 
els. This combines the high accuracy on large scales ('^ 1%) 
provided by systematic perturbative expansions with the 
reasonably good accuracy on small scales (~ 10%) pro- 
vided by the halo model. This is particularly useful on the 
weakly nonlinear scales associated with the baryon acous- 
tic oscillations, where the matter power spectrum and bis- 
pectrum show small "wiggles" that would be difficult to 
reproduce with a high accuracy by a simple halo model 
or fits to simulations (unless one runs a simulation with 
the required cosmological parameters). As a first step in 
th is direction, we have recently built such a unified model 
in IValageas fc Nishimichil ()2011h . focusing on the matter 
power spectrum. In this paper, we extend this study to the 
matter bispectrum, using the same halo model and pertur- 
bative schemes. As for the power spectrum, we explain that 
the "1-halo" and "2-halo" contributions contain new coun- 
terterms that ensure a physical behavior on large scales. 
Then, we show that we obtain a good agreement with nu- 
merical simulations without introducing new parameters. 
We also note that it is possible to improve the predictions 
for the power spectrum on mildly nonlinear scales using the 
shape of the predicted bispectrum. 

This paper is organized as follows. We first describe in 
Sect. [5] the computation of the density bispectrum from a 
Lagrangian point of view, making the link with halo mod- 
els and perturbative expansions through the decomposi- 
tion into 1-halo, 2-halo, and 3-halo contributions. Then, 
we present in Sect. |3] a simple implementation, based on 
the halo model and the pert urbative approaches used in 
IValageas fc Nishimichil ()201lD for the power spectrum. We 
compare our results for the bispectrum with numerical sim- 
ulations in Sect. SI for equilateral and isosceles configura- 
tions, from linear to highly nonlinear scales. Next, we dis- 
cuss the differences between our work and some previous 
studies in Sect. jSJ and we investigate in Sect. [Hj the depen- 
dence of the predictions on the choice of the perturbative 
scheme and on the halo parameters. We explain in Sect. [7] 
how to use the shape of the predicted reduced bispectrum 
to improve the predictions for the power spectrum and the 
bispectrum, and we estimate the accuracy of our unified 
model in Sect. [5] before concluding in Sect. [HI 



2. Decomposition of the density bispectrum from a 
Lagrangian point of view 

We show in this section how the density bispectrum can 
be split into 3-halo, 2-halo, and 1-halo terms from a 
Lagrangian point of view, and we make the connection with 
perturbation theory. 

2.1. Lagrangian-space formulation 

In a Lagrangian framework, one considers the trajectories 
x(q, t) of the particles, of initial Lagrangian coordinates q, 
and Eulerian coordinates x at time t. At any time t, this 
defines a mapping, q i— > x, from Lagrangian to Eulerian 
space, which fully determines the Eulerian density field p(x) 
through the conservation of matter, 



p(x) dx = pdq. 



(1) 



where p is the mean comoving matter density of the 
Universe and we work in comoving coordinates. Then, 
defining the density contrast as 



P(x,i) - p 



P 



and its Fourier transform as 
dx 



'^"^ J (2.)3 
we obtain from Eq.([T]) 



^(x). 



<5(k) 



dq 

(27r)3 



g-ik x(q) _ g-ik q 



(2) 



(3) 



(4) 



We define the power spectrum, P{k), and the bispectrum, 
B{ki,k2,k3), which are the Fourier transforms of the two- 
point and three-point density correlation functions, by 

(5(ki)5(k2)) = (5D(ki+k2)F(fci), (5) 

(5(ki)5(k2)5(k3)) = Soihi + k2 + kg) S(fci, k2, k^) . (6) 

Here we used the fact that the system is statistically homo- 
geneous, which gives rise to the Dirac prefactors associated 
with statistical invariance through translations, and statis- 
tically isotropic, which implies that P{ki) and B{ki, ^2, fcg) 
only depend on the lengths |ki| and {|ki|, |k2|, Ikaj}. 

We have already studied the power sp ectrum in a pre- 
vious work (j Valageas fc Nishimichil [201 it) , so that we here 
focus on the bispectrum. Our goal is to follow the same 
approach to combine the perturbation theory and the halo 
model and to build a simple unified model. Within this 
framework, we wish to decompose the bispectrum into "per- 
turbative" and "non-perturbative" contributions. The for- 
mer is associated with configurations of particle triplets 
(or p— uplets for the p— point correlation function) where 
the particles belong to different halos, whereas the lat- 
ter is associated with configuratio ns where several parti- 
cles b elong to the same halo. As in IValageas fc Nishimichil 
(|201lD . our strategy is to evaluate the perturbative contri- 
bution through perturbation theory (either with the stan- 
dard expansion or with resummation schemes) and the non- 
perturbative contribution through a halo model. Then, the 
Lagrangian framework allows us to recover the countert- 
erms that arise in the halo-model contributions and that 
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are necessary to ensure a good behavior on large scales. In 
addition, it is of interest to see how the standard expression 
for the bispectrum obta ined from an Eulerian point of view 
(iCoorav fc Shethl[200l can be recovered (with the addi- 
tion of these new counterterms) from a Lagrangian point 

of view^ 

In IValageas fc Nishimichil ()2011[ ). where we studied 
the power spectrum, we could factorize out the Dirac 
prefactor of Eq.© and write for the power spectrurn 
the well-known expressio n ([Schneider fc Bartelmannlll995l: 
iTavlor fc Hamiltonlll996[l 



P{k) = 



dq 

(27r)3 



^-ik-Aj 



pikq 



(7) 



where we introduced the Eulerian-space separation Ax — 
x(q) — x(0). Thus, Eq.([7]) only depends on the relative dis- 
placements of the particles, as it must because uniform 
translations do not affect the structural properties of the 
density field. Then, Eq.© provides a convenient starting 
point because the Dirac prefactor of Eq.© has already 
been taken into account and it will not be put in danger by 
approximations used in later steps. 

For the bispectrum, it is still possible to factorize out 
the Dirac prefactor of Eq.® by introducing the relative 
displacements. This yields 



B{ki,k2,kz) 



dq2q3 

(27r)6 



) 



-,5i3(ki)F(fc2) - 5D{^2)P{ki) - (5D(k3)P(fcl) 

-,5i3(ki)(5o(k2), (8) 

where we introduced Ax2 = x(q2) — x(0) and Ax3 — 
^(qa) ~ ^(0)- However, the symmetry over {fci,fc2,fc3} is 
no longer transparent in Eq.([5]) and it is not very conve- 
nient to count the triplets that are within one, two, or three 
halos. Therefore, we prefer to work directly with the third- 
order moment although this will require an additional 
approxi mation to those needed t o com pute the power spec- 
trum in IValageas fc Nishimichil ()2011[) . Thus, using Eq.(|3]) 
we write 

dqidq2q3 



((5(ki)(5(k2)(5(k3)) = 



(27r)s 



X e 



-ik. 



where Xj — x(qj). Next, we split the average ^ into three 
contributions, associated with the cases where the triplet 
{qijq2,q3} belongs to either one, two, or three halos: 

(^(ki)^(k2)5(k3)) = (J(ki)5(k2)5(k3))iH 

+ (^(ki)^(k2)J(k3)}2H + (^(ki)5(k2)5(k3))3H. (10) 

This assumes that all the matter is contained within halos, 
so that these three contributions sum up to Eq. ^ , counting 
all particles. 

2.2. "1-halo" contribution 

Let us first consider the 1-halo contribution. It can be writ- 
ten as 

3 

(^(ki)J(k2)<5(k3))iH = (/n (e-ikrx. _e-ik,.q,) 

3 

a j=l 



(11) 



where the sum runs over all halos, labeled by the index a, 
and the top-hat factors 6{qj G Ma) are unity if the particle 
qj belongs to the halo a and zero otherwise. This clearly 
gives the contribution to Eq.® of the configurations where 
the triplet {qi,q2,q3} belongs to a single halo. Equation 
([TT|) also reads as 



(J(ki)^(k2)^(k3)) 



n 



dq 



J_ /'p-ikj Xj _ -ikj-qj 



11 (27r)3 V 

3 

J dq^dA^n(q^M) Heiq^eM)), 



(12) 



where n(q^, M) is the halo number density in a given real- 
ization of the initial conditions (as denoted by the hat, in 
contrast with its mean). 



n(q^ M) = J2 Soicf - q^J Sd{M - M„). 



(13) 



Here q^ is the "center" of the halo a, which may be de- 
fined as the halo center of mass, in Lagrangian space (i.e., 
in the initial or linear density field). Note that in Eg. dT^ 
we have chosen to count halos in Lagrangian space, by their 
Lagrangian center q^, but we could as well count them in 
Eulerian space by using the Eulerian-space center of mass 
x^. The mean of the halo number density does not de- 
pend on location, thanks to statistical homogeneity, and it 
is given by the halo mass function. 



{n{q',M))^n{M), 
which we write as 

r.(M)dM = ^/(^)-, with ^ = -^. 



(14) 



(15) 



As usual, we have introduced in Ea. lfT^ the scaling function 
j(y) and the reduced variable z^, where ct(M) is the rms 
linear density contrast at scale M, or Lagrangian radius 
qm, with 



47r 

ct(A/) = (j{qM) with M = p—QM, 
and 



a^iQM) = 4^ / dfcfc2Pi(fc)VF(fcgM)' 



(16) 



(17) 



where W{kqM) is the Fourier transform of the top-hat of 
radius qm, defined as 

WikqM)^[ ^e'--^ ^ 3 Bin(fcgM)-fegM C0s(fcgM) _(^3) 
Jvm Vm (kqMr 

In the second Ea. (|T5|) . the linear density contrast Sl is 
related to the nonlinear density threshold S that defines 
the halos t hrough the sph erical collapse dynamics, as (5 = 
T{5l), see lValageasI (|2009D . Thus, Eq.^ also writes as 

(^(ki)^(k2)^(k3))iH = /dq^VM^^"^^ 



)) 



qSAf, 



(19) 
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where the average (..)qc^M is the conditional average un- 
der the constraint that the three particles belong to the 
same halo of center q'', mass M, and L agrangian volume 
Vm- As in lValageas &: Nishimichil ()2011[ ). we make the ap- 
proximation of fully virialized halos, that is, the particle q 
has lost all memory of its initial location and velocity and 
it is located at random within the halo. This means that in 
Eulerian space the average is defined by the density profile 
of the halo, Px''.m{^), as 



.M 



1 

M 



dxe-''-''pAf(|x-x^|), 



(20) 



because q^ can be identified with a uniform probability to 
all the particles that make up the halo (approximation of 
complete relaxation). Here we also made the approxima- 
tion of spherically symmetric halos, and introducing the 
normalized Fourier transform of the halo radial profile. 



UM{k) 



^-Ta =M^^^ Pm{x), 



we obtain 



(e-'''-''^)qSA/ = e-"'-^^ UM{k). 



(21) 



(22) 



Making the approximation that these halos are also spher 
ically symmetric in Lagrangian space, we can use Eq.([T 
for the factors e~'^^^''^^ and Ea. (fTO|) reads as 

(J(ki)^(k2)^(k3))iH = / dq^— ^/(;.)e-('^^+''^+'^3) q' 



V AT 



M 



p(27r)2 



3 3 



n (e-*^-*^{iAf(fc,) - W{k,qM)) , (23) 



where we introduced the displacement, = x*^ — q'^, of 
the center of mass of the halo. We have written Eas. (PD)) 
and ((22|) - (|23)) as if the relation x'^(q'^) were determinis- 
tic, but a priori we should take the average over the dis- 
placement in Ea. (P5)) . First, we note that (e"'''^ *'') 
does not depend on q'^, thanks to statistical homogene- 
ity, so that the integral over q*^ in Eg. ([25)1 gives the ex- 
pected Dirac factor (5£i(ki -I- k2 -I- ks), in agreement with 
Eq.([6]). Then, the displacement only appears in mixed 
terms, such as e~'^'^i+'^^^'*''uM(^i)'iiAf(fe)W^(fc3'?Af)- Going 
back to Eq.([3]), we can see that the term W{k^qM) arises 
from a factor (ka), so that in the regime where the 1- 
halo contribution is dominant, that is ((5(ki)(5(k2)i5(k3)) ~ 
((5(ki)(5(k2)(5(k3))iH, it should reduce to SdO^s). In Eq.iP]) 
this corresponds to the fact that \W{k^qM)\ ^ 1 for k^ :s> 
1/qM- Then, thanks to the prefactor (5D(ki-|-k2-l-k3) we can 
see that in this regime |ki -f k2| — ?• and i 
In order to satisfy these properties, we simply make the ap- 
proximation = 0, that is, we neglect the displacements 
of halos. Then, performing the integral over q'^, which al- 
lows the factorization of the Dirac factor S^Cki + k2 -I- k3), 
we obtain the 1-halo contribution to the bispectrum as 



di^ 



BlH(fcl,fc2,fc3) = / — /H 1 ^(27r)3 



M 



V 

3 

X J]^ (uMikj) - W{kjqM 



(24) 



Thus, as for th e pow er spectrum studied in 
IValageas fc Nishimichil (|201l[ ). we recover the stan- 
dard expression of the 1-halo contribution to the 
bisp ectrum derived within an Eulerian fra m ework , 
see ICoorav fc ShethI (|2002D : IScoccimarro et al.l (|200l[ ). 
with the addition of the new counterterms associated 
with the factors W . No te that for the pow er spec- 
trum we obtained in iValaeeas fc Nishimichil (|201ll) 
a factor {uM{k)^ — W{kqM)'^^, instead of the factor 
[MM(fc) — W^(fc<ZAf)]^ that we would obtain using the 
approach described above. The difference arises because 
for the power spectrum we used Eq.©, where the Dirac 
prefactor bu^\ +k2) has already been factorized, whereas 
for the bispectrum we used Eq.® instead of Eq.®. The 
advantage of Eqs.© and (O is that they only involve 
relative positions. Ax, so that Eulerian position x'^ of the 
halo never appears (for particles located in the same halo) 
and we do not need to introduce any approximation for 
the halo displacement . 

However, for higher-order functions, such as the bis- 
pectrum or the trispectrum, the symmetric equation ^ 
(and its A^-point generalization) provides a simpler start- 
ing point. In particular, the counting of the volumes over 
{qi, .., qTv} associated with a single halo simply factorizes 
as y^ii whereas introducing relative positions makes the 
geometrical countings slightly more intricate. 

The factors [uAf(fcj) — W^(fciQA/)] cannot be interpreted 
as modifications of the halo profiles puix) under the form 
of compensated profiles (the result pu [x) — p would actu- 
ally be negative at large x). Indeed, both quantities do not 
have the same radius in real space, and actually apply to 
two different spaces, namely the Eulerian and Lagrangian 
spaces. 

A second important feature is that the expression (1111) 
automatically ensures that the 1-halo contribution decays 
at least as i3iH(fci, ^2, fcs) (x fcj for any kj 0, since no 
linear dependence on kj can remain because of statistical 
homogeneity and isotropy. This can be checked on the final 
expression (|24t . since ua^(O) = M^(0) = 1 and the func- 
tions UA/(fc) and Wikqu) are regular functions of fc^ at the 
origin. 



fcj -> : 



SiH(fci,fc2,fc3) fc|. 



(25) 



This agrees with the requirements associated with a small- 
scale redistribution of matter. Here, there is a simplification 
as compared wi th the m a tter p ower spectrum P(fc). Indeed, 
as discussed in iPeeblesI (|1974[ ). small-scale redistributions 
of matter generically give a low-fc tail ^(k) oc fc, whence 
P(fc) oc fc^, while taking into account momentum conserva- 
tion implies the steeper d ecay J(k) oc fc'^ , whence -P(fc) oc 
fc . Thus, in Valag eas & Nishimichil (poll we recovered the 

expected fc^ tail, through the factor (uA/(fc)^ — W^(fc<ZA/)^), 
since our analysis did not enforce momentum conservation. 
Using the approach described above, we would actually ob- 
tain a fc'' tail, through the factor [MA^(fc) — W^(fc(ZA/)]^, once 
we make the approximation '^'^ = 0. This clearly satis- 
fies momentum conservation (the peculiar momentum of 
each halo is zero) but comes at the expense of an addi- 
tional approximation and one should not give too much 
weight to this property. By contrast, for the bispectrum 
(and higher-order correlations) small-scale redistributions 
of matter lead to the same fc| tail, whether we take into 
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account momentum conservation or not, because the lin- 
ear term kj vanishes by symmetry. Note that when all 
wavenumbers go to zero at the same rate, we have from 
Ea. ([M)) the decrease 



A^O: BiH(Afci,Afc2,Afc3) oc 



(26) 



As noticed in previous works (jCoorav fc ShethI 1200 2V in 
the usual version of the halo model the 1-halo contribu- 
tion to the bispectrum does not contain the counterterms 
W{kjqM), so that it goes to a strictly positive constant 
at large scales and dominates over the lowest-order pre- 
diction of perturbation theory. This led to rather inaccu- 
rate predictions on large scales (an overestimate of 20% at 
k O.Ol/iMpc^^), because for CDM cosmologies the per- 
turbative prediction scales as Sport ~ PL{k)^ oc with 
Us — 1 (for kj ~ k for all j). As seen from Eas. (P5)) - (P5)) . 
our approach solves this problem (we will return to this 
point in Sect. [5]). 



2.3. "2-halo" contribution 

We now turn to the 2-halo contribution to the bispectrum, 
following the procedure described in the previous section 
for the 1-halo term. Thus, as in Ea. lfT^ we can write the 
2-halo contribution as 



3 

(^(ki)^(k2)<5(k3))2H = (y n 



L (2^)3 



— ik^ 



X / dqidMidq2dM2n(qi,Mi)n(q2,M2) 

X0(qi e Ah) e{(i2 e Ah) 0(q3 G Ah) ) + 2 eye, (27) 

where we sum over all pairs of distinct halos, of mass Mi 
and M2, which contain one and two of the particles q^ . Here 
we have written the contribution where the halo AIi con- 
tains the particle qi, and the two additional contributions, 
noted "2 eye." , correspond to the cases where AIi contains 
the particle q2 or qa. Using again the approximation (22) 
of fully relaxed halos, the average of Ea. (E71) writes as 

(<5(ki)^(k2)^(k3))2H= / dq5^dq^2— T^fi'^i)fi^2) 



Ah Mi 



vi V2 A/hA'h 

iki-q? -i(k2+k3)-q2 



X (e-*^-*i«Af,(fci)-W^(fci(ZMi)) 
3 

X n (e""'^'*^ UM^kj) ~ W{k,qM2)) + 2cyc., (28) 
i=2 

where fi2(|q2 ~ Qil) is the two-point correlation function 
of halos of mass Ah and Ah, in Lagrangian space. As a 
first step, let us neglect halo motions, '4''^ = 0. Then, as 
expected, the contribution associated with the factor 1 in 
(1+^12) vanishes (because the integral over q^; yields 6d{}^i) 
and [uMi(fci) — W{kiqi,[i)] = for ki = 0), so that only 
the term associated with the large-scale correlation of halos 
remains. Writing this two-point correlation in terms of the 
halo power spectrum. 



the integrals over q^ and give the Dirac factors 5d (k + 
ki)(5£)(k — k2 — ka). Therefore, we recover the Dirac pref- 
actor ^^(ki -|- k2 -I- ka), which we can factorize out to write 
the 2-halo contribution to the bispectrum as 

S2H(fcl,fc2,fc3) = / — — ^^I{v,)f{v2)Pl2{kl) 

3 

X (uMiiki) - W{kiqMi)'j Y[ {^Miikj) ~ W{kjqM2)^ 
+2 eye. (30) 



Ci2(|q§ - q^l) - J dke^'-^'i^-i?) P^k), 



(29) 



Then, as is customary, we could write the halo power spec- 
trum as Pi2(fc) = &(Mi)&(M2)P(fc), where b{M) is a mass- 
dependent bias factor (of order unity for typical halos) that 
we approximate as scale-independent (thereby neglecting 
exclusion constraints between halos). 

However, the expression (15(1)) is not really satisfactory. 
Indeed, for ki it decreases as kfP{ki) because of 
the prefactor [uMiiki) — W{kiqMi)], while we would ex- 
pect to recover a behavior proportional to P(fci). Indeed, 
in this regime this 2-halo contribution should describe the 
power that arises from the large-scale correlation between 
density fluctuations at a position qi and at a position 
q2, with |q2 — qi| ^ ^ir/ki, the fluctuations at q2 be- 
ing furthermore decomposed over smaller-scale fluctuations 
within single halos of mass M2. In other words, the large- 
scale power P{ki) should not be damped by the pref- 
actor [uMi(fci) — W{kiqMi)]- Going back to Ea. (l28)) . we 
can see that the product w^/i (fci)wM2 (fc2)MM2 (^3) involves 
a prefactor e'''^'*-*^-*!)^ which only depends on relative 
displacements and as such is physically meaningful and 
should be taken into account. This should be contrasted 
with the behavior obtained for the 1-halo contribution ((23|) . 
where the product of three terms um gave the prefactor 
gi(ki-i-k2+k3) * ^ -^vjiich had to disappear because it does not 
depend on relative displacements, and was indeed irrelevant 
thanks to the Dirac prefactor (5_D(ki -|-k2 -l-k3). This means 
that the approximation '5''^ = is not as good for the 2- 
halo contribution as for the 1-halo contribution, because it 
neglects relative halo displacements, which did not appear 
in the latter case (i.e. its consequences are stronger in the 
former case). 

In order to have a 2-halo contribution that behaves in a 
reasonable fashion, we choose to use instead of Eg. ([50)) the 
simple expression 

P2H(fcl,fc2,fc3) = PL(fcl) / — -^/(i^) 

J V p(27r)-^ 

3 

X n ("m(%) - WikjqM)) + 2 eye. (31) 

Thus, we have made the approximation Pi2{k) ~ PL{k), 
where Pbik) is the linear matter density power spec- 
trum, and we have set the spurious prefactor [{(Mi(^i) — 
W{kiqMi)] to unity. In principle, instead of the linear power 
Phik) we could include higher orders of perturbation the- 
ory. However, in view of the approximate nature of Ea. (l3ip 
this is not really necessary, especially since the 2-halo con- 
tribution is subdominant on most scales of interest, as 
shown by the numerical results in Sect. |3| 
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In the large-scale limit we obtain from Eg. ipTj) the 
asymptotic behaviors 

fe, ^ : B2H(fci, k2,h) (X PL{k,), (32) 

when only one wavenumber goes to zero, and 

A^O: S2H(Afci,Afc2,Afc3) oc A^Fl(A), (33) 

when all wavenumbers decrease at the same rate. In Ea.(|5^ 
we used the fact that for CDM cosmologies we have Pl {k) ~ 
/c"^ with ris ~ 1 at low fc. It is interesting to note that the 
behavior p2p agrees with perturbation theory, as shown 
in App. El and Ea. (|A.5p . There, we note that this scaling, 
which applies when only one of the wavenumbers kj goes 
to zero, for instance ki , is valid at all orders of perturba- 
tion theory. In this regime, the 2-halo contribution, which is 
non-perturbative (as discussed in Sect. 12.41 below), is there- 
fore on the same order (i.e., ^ Pl(A;i)) as the perturbative 
contribution. This can be understood from the arguments 
developed above to obtain the expression (PT|) . as the non- 
perturbative effects associated with the formation of small- 
scale nonlinear structures around the nearby points X2 and 
X3 (or q2 and qa in Lagrangian space) do not modify the 
larger-scale correlation with the very distant point Xi . That 
is, the non-perturbative effects associated with the forma- 
tion of a halo correspond to a local redistribution of matter. 
Then, at leading order, non-perturbative terms simply lead 
to a "renormalization" of the coefficient that multiplies the 
factor Phiki), by an amount that depends on how far in 
the nonlinear regime the wavenumbers ^2 and fcs are. 

When all wavenumbers go to zero, the property 
again ensures that the 2-halo contribution becomes negligi- 
ble as compared with the lowest order term of perturbation 
theory, which scales as Pl{X)^. Thus, as for the 1-halo con- 
tribution, the counterterms W of Ea. ((3T|) solve the prob- 
lem encountered with the usual implementation of the halo 
model, which does not recover perturbation theory on very 
large scales. 

2.4. "3-halo" contribution as a perturbative contribution 

For the 3-halo contr ibution we follow the sp i rit of the ap- 
proach described in IValaeeas fc Nishimichil (|201lf ) , as we 
bypass the halo model to make contact with standard 
perturbation theory. Indeed, as we have seen above for 
the 2-halo term, it is not straightforward to take into ac- 
count large-scale correlations within this Lagrangian im- 
plementation of the halo model, starting with the expres- 
sion (|9]). This would require some modeling of relative dis- 
placements, which may not be worth the effort (if one 
wishes to improve in a significant manner over the ap- 
proximation (j34p used below). Thus, we take a much more 
direct and simpler route, noting that the 1-halo and 2- 
halo contributions vanish at all orders of perturbation the- 
ory. Here, we mean by perturbation theory any expan- 
sion over powers of the linear power spectrum P^ (for 
the Gaussian initial conditions that we consider in this 
article), such as the standard perturbation theory derived 
within t he fluid approximatio n to the equati ons of motion 
(|Goroff et al. 1986; Bcrnarde au et all l2002ah . Indeed, the 
halo mass function shows a large-mass tail of the form 
e~'^ = e'^^L^^'^" (^^)), which vanishes exponentially for 
Pl ^ and has a Taylor expansion over powers of Pl that 




Fig. 1. Probability F3H that the triplets associated with 
the equilateral bispectrum Beq{k) belong to three different 
halos, from Ea.([55|). at redshifts z — 0.35, 1, and 3. 

is identically zero. Therefore, the remaining 3-halo contri- 
bution to Eq. (|T0| must be consistent with perturbation the- 
ory up to all orders, and we simply use the approximation 

Bmikl,k2,k3) = Sport (fci, ^2, ^3), (34) 

where Bpert is the bispectrum obtained from perturba- 
tion theory. To follow more closely the approach used in 
IValageas fc Nishimichil (|2011[ ) for the power spectrum, we 
should multiply iJpcrt in Ea. ((M)) by a factor F^h, with 
< < 1, which counts the fraction of triplets that 
are located within three distinct halos. From the same ar- 
guments we have F^n = 1 at all orders of perturbation 
theory, and it only differs from unity by non-perturbative 
terms such as e~*J-/'^°' where M is typically the mass 
scale associated with the maximum of {27T/kj}. For illus- 
tration purposes, we display in Fig.[T]the following estimate 

of i^3H, 

FMk)=(^£^f{i^?l\ (35) 

where Vk = ^L/f^iqk) with = 27r/fc. From Eq. ljlSp . this is 
the probability that three particles belong to halos of size 
smaller than g^, neglecting halo correlations, and it should 
give an estimate of the range where i^an — 1 (for the equi- 
lateral bispectrum, which involves a single wavenumber k). 
As expected, the comparison with Fig. [5] below shows that 
the departure from unity of F3H roughly corresponds to the 
scale where the 1-halo and 2-halo terms become of the same 
order as the 3-halo term. In view of the approximations in- 
volved in these 1-halo and 2-halo contributions, we do not 
try to include the deviations from unity of -Fsh, which only 
play a role in the transition regime. Therefore, we make 
the approximation i'sH — 1 and use the simple prescription 
(El- 

In practice, the perturbative term Bpert is not exactly 
known (assuming the perturbative series is convergent), 
and one must truncate the perturbative expansion. There 
are many ways to do so, because a priori one can use 
the standard perturbation theory or any alternative expan- 
sion scheme, which corresponds to partial resummations of 
the st andard diagrams. As shown in I Valageas fc Nishimichil 
(|201lD . for the power spectrum it happens that standard 
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perturbation theory is not a good choice, because higher or- 
der terms grow increasingly fast on small scales and prevent 
a good description of the highly nonlinear regime. Then, 
one must use resummation schemes that agree with stan- 
dard perturbation theory up to the required order while 
remaining well-behaved on small scales. As we will check 
in Sect. SI for the bispectrum it turns out that at one-loop 
order the standard perturbation theory remains small at 
high k compared with the 1-halo and 2-halo contributions, 
which makes this an efficient choice; but we will also study 
the direct steepest-descent resummation scheme described 
in detail in I Valageas fc Nishimichil (j201lh for the computa- 
tion of the power spectrum. 

In most studies that are based o n the halo model , 
the 3-halo term is rat her written as dMa fc Frvl [2"000bt 
IScoccimarro et al.ll200ll : ICoorav fc Shethll2002f ) 



l[{b{M)aMih)) 



Strco(fcl,fc2,fc3), (36) 



where {b{M)uM) is the average over mass, weighted by the 
halo mass function, of the bias parameter b and halo profile 
um] -Btree IS the matter bispectrum obtained at lowest or- 
der through perturbation theory, given by Eq. (j39p below. In 
our approach, we do not introduce this halo bias parameter 
because we simply write the 3-halo term as the perturbative 
matter bispectrum (IMl) . This implies that we focus on the 
matter bispectrum because our approach does not contain 
the halo bispectrum B^^^- {ki^k^jk^; Mi, M3) as an in- 
termediate tool. Therefore, to compute the bispectrum of 
halos of given masses we would need to add another pre- 
scription. On the other hand, this makes our model for the 
matter bispectrum simpler and more robust because it does 
not rely on bias parameters that are fairly difficult to com- 
pute in a systematic and well-controlled manner (because 
halos themselves do not appear in a natural fashion in the 
equations of motion). Moreover, Eg. dMl) ensures that our 
results agree with perturbation theory up to the order of 
truncation. In this paper, we will go up to 1-loop order, 
while most previous studies that involved the halo model 
only used the tree-level contribution, as in Eq. (|36t . 

2.5. Comparison with the Eulerian framework 

As seen in the previous sections, the derivation of the bis- 
pectrum from a Lagrangian implementation of the halo 
model is not as straightforward nor as "clean " as for the 
power spectrum (| Valageas fc Nishimichil [20TTh . Indeed, we 
had to introduce some additional approximation regarding 
the halo displacements and for the 2-halo term that 
mixes larger-scale and intra-halo wavenumbers we had to 
partially bypass the naive prediction of this halo model to 
recover the large-scale behavior ([5^ . 

More generally, it is not surprising that the Lagrangian 
implementation of the halo model requires more steps than 
the usual Eulerian version. Indeed, within an Eulerian 
framework we only need to describe the density field at 
the time of interest t. Within a Lagrangian framework, we 
need to add some information on the building of this density 
field, that is, the mapping between the initial coordinates 
q and the current positions x of the particles. This neces- 
sarily implies that a practical Lagrangian implementation 
requires more hypothesis or approximations (because there 
are more intermediate quantities, even though in principle 



the results would be the same if we made no approximation 
at all). However, we think it remains of interest to describe 
how TV-point correlation functions or poly-spectra can be 
constructed within such a Lagrangian framework, based on 
the halo model, as we have done above. First, it is use- 
ful to see how identical or similar approximate models can 
be derived from different points of view, since this gives 
more weight to the results. Second, it provides additional 
insight on the underlying approximations and it could of- 
fer an alternative route to more precise modeling. Third, 
it provides a natural derivation of the counterterms W in 
the 1-halo contribution ([M)) . which ensure the large-scale 
behaviors (P5)l - (|26p that were missed in previous Eulerian 
implementations. 

As explained above, for the 2-halo contribution we had 
to take some liberty with a naive implementation of the halo 
model to recover the asymptotic scaling (15^ . However, in 
view of the approximate nature of these Lagrangian halo 
models, we think it is best to be guided by physical argu- 
ments and to make sure physical constraints are satisfied, 
instead of strictly adhering to approximate models that can 
show a varying degree of accuracy, depending on the quan- 
tities under scrutiny^ 

However, as in Valageas fc Nishimichil (|201l[ ). it is in- 
teresting to note that the counterterms of Egs. fM)) and 
([5T|) might also be guessed within an Eulerian framework, 
from the requirement that the bispectrum should vanish for 
a perfectly homogeneous universe. Indeed, with the usual 
halo-model approximations, we can still split a constant- 
density medium over the same set of "halos", or cells, 
of radius qM and center q'^ (neglecting geometrical con- 
straints associated with the impossibility of covering a 3D 
space with spherical cells). The only difi'erence is that be- 
cause these "halos" have not changed and have remained 
at the constant density p, their Eulerian and Lagrangian 
properties are the same. In particular, the Eulerian radius 
xa/ is equal to qm, and the halo profiles are simply given 
by UM{k) = W{kqM) for this uniform system. Then, the 
counterterms in Eas. ([Ml) and (PT|) clearly satisfy the con- 
straint _B(fci,fc2,fc3) = for this homogeneous universe. 
However, by itself this argument is not sufficient to imply 
the precise form of Eg. ([Ml) (nor of Eq . ([?!]) ') . since the choice 

wm(%) — rij W{kjqM) would also satisfy this constraint 
(but not the properties ([^51 ). 

Although the main reason for taking into account the 
counterterms W in the 1-halo and 2-halo contributions 
is to satisfy physical requirements, this is also needed 
to reach hig h accu racy on weakly nonlinear scales. Thus, 
iGuo fc Jind ()2009f ) find that on scales on the order of 
0.1/iMpc~^ their 1-halo and 2-halo terms (which do not 
include these counterterms) are still too large and spoil the 
agreement of (tree-order) perturbation theory with numer- 
ical simulations. Giving a faster decay on large scales of 
these 1-halo and 2-halo contributions, the counterterms W 
improve the agreement of the analytical predictions with 
the simulations, as we will check in Sect. [S] 

3. A simple implementation 

We briefly describe in this section the numerical implemen- 
tation of our model for the matt er density bispectrum. 
Furth er details can be found in IValageas fc Nishimichil 
(|20T1 . 
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3.1. Halo properties 



For numerical comp utations we u se the sa me ha lo model 
as the one used in iValaeeas fc Nis himichi (|201ll ) for the 
ower spectrum . Thus, the halo mass function is given by 
Valageasl[20n9h 



i 



f{v) = 0.502 [(0.6 lyf-^ + (0.62 z.)"-^] er"'^^, 



(37) 



which is normalized to unity and provides a good match to 
numerical simulations. We choo se the usual NFW pr ofile for 
the halo density profile pm{x) (jNavarro et al.llT997h . which 
also has an expli cit form for its Fourier transform UM{k) 
(jScoccimarro et a l. 2001). In particular, we truncate the 
halo profiles at the density contrast 5{< r2oo) = 200, which 
defines their radius r2oo. For the concentration parameter 
we take 



C(M2( 



10.04 



M2 



2 X 10l2/i-lM(r 



(1 + ^) 



-0.8 



(38) 



which is si milar to the behaviors measured in numerical 
simulations (iDolag et all 120041: iDufltv et aT]|2008D . and was 
found in IValagcas fc Nishimichi ( 201lt) to provide a good 
tool for the density power spectrum (in this sense, it would 
also describe to some degree the effect of halo substruc- 
tures). 

3.2. Perturbative contribution 

For the perturbative contribution ([M)) to the bispectrum, 
we investigate both the standard perturbation theory and 
the "direct steepest-descent" resummation scheme intro- 
duced (among other perturbative expansions) in iValageas 
(jlOOT^)- We only go up to 1-loop order for both methods, 
so that the first approach only contains all 1-loop diagrams, 
while the second approach also includes partial resumma- 
tions of diagrams at all orders. 

A detailed description of these m ethods fo r two-p oint 
and three-point functions is given in IValageasI (|2008l) . To 
facilitate the theoretical comparison between both ap- 
proaches it is convenient to write them with the same di- 
agrammatic language. Then, the diagrams associated with 
standard perturbation theory are shown in Fig. [5J where 
the solid lines are the linear two-point functions, either the 
linear correlation (lines without arrows) or the linear re- 
sponse (lines with an arrow that shows the fiow of time). 
The three-leg vertex is the kernel Kg that appears in the 
equation of motion, which can be written as O-V' = Ka-^^ip^ 
where C is a linear operator and t/i is a two-component vec- 
tor that describes both the density and velocity fields. The 
diagrams shown in Fig. [5] are not those associated with the 
usual description of standard perturbation theory, which 
is described in App. \^ but of course they are equivalent. 
The expansion of Fig. [2] applies to the three-point correla- 
tion C3 = (ipipip), which contains the density and velocity 
three-point functions, as well as their cross correlations, but 
in this article we only consider the density three-point cor- 
relation, that is, the matter density bispectrum in Fourier 
space. 

The usual description of standard perturbation theory 
arises from the expansion of the density and velocity fields 
over powers of the linear field Sl, as in Eqs. (IA.l[) - (|A.2p . 
Then, the TV-point correlations are obtained by taking the 
Gaussian average of the product of these expansions, as in 



+ 48^ 



+ 24 



+ 48 



+ 48 




+ 24 — ^ 




"standard" 



Fig. 2. Diagrams associated with the standard perturba- 
tion theory as obtained from a path-integral formalism up 
to 1-loop order for the three-point correlation C3 . Although 
they are different from those obtained by the standard ap- 
proach, their sum at each order is identical to the sum of 
the standard diagrams of the same order over P£,(fc). The 
lines are the linear two-point correlation Cl (blue solid line) 
and the linear response Rl (red solid line with an arrow). 
The black dots are the three-leg vertex Kg that enters the 
quadratic term of the equation of motion. The numbers are 
the multiplicity factors of each diagram. The tree-order di- 
agram (a) gives Ea. ([M)) . for the density bispectrum, while 
the 1-loop diagrams (b),..,(i), give the contribution of order 

Pi 



Ea. (|A.3[) . using Wick's theorem as in Ea. (|A.4[) . This yields 
diagrams with vertices Fn that have an increasing number 
n of legs as one goes to higher orders, see also Fig. lA.ll 
The diagrams of Fig. [5] are obtained from a path-integral 
formalism, where the averages over the initial conditions 
are already taken and one directly works with correlation 
and response functions. Then, expanding over powers of 
the nonlinear interaction vertex Kg gives the expansion of 
Fig. [21 It looks quite different from the standard diagrams, 
since only one vertex appears, i.e., the three- leg kernel i^T^, 
whatever the order of the expansion. On the other hand, 
in addition to the linear power spectrum these diagrams 
also involve the linear response function. Then, the order 
of the expansion corresponds to the number of loops of the 
diagrams. 

Since both expansions can also be written as expansions 
over powers of the linear power spectrum Pi(fe), they are 
actually equivalent, at each order. However, at a given or- 
der, there can be a different number of diagrams between 
both expansions, and it is only the two sums over all dia- 
grams of that order over Phik) that coincide. In particular, 
it can be seen that diagrams (h) and (i) in Fig. [5] show 
an ultraviolet divergence for linear power spectra with a 
large-fc tail Pl(A;) oc fc" with n > —3. However, these two 
divergences cancel out and the sum is finite for n < — 1, as 
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in the case of the standard diagrams (see IValageasI (|2008[ ) 
for details). 

The lowest-order contribution, or "tree-order" term, is 
given by diagram (a), which yields the well-known result 



, 4(k2 • ks f 



"10 


(k2 h\ 







+ 2 eye. 



^2 ■ k3 

(39) 



The explicit expressions of the 1-loop order diagrams 
(b),..,(i), which are on the order of P| (as shown by the 
three bl ue solid lines w hich they contain), are given in 
App. A of lValageasI EoM) . For numerical computations it is 
convenient to perform analytically integrations over angles 
instead of directly implementing these expressions into the 
codes. Although this requires some care (since angular in- 
tegrations may yield trigonometric or hyperbolic functions, 
depending on the amplitude of the wavenumbers) , there is 
no fundamental difficulty. 



Ca= 6 




+ 



+ 48 



+ 48 




+ 24 




"direct 
steepest-descent 
resummation" 



Fig. 3. Diagrams obtained at 1-loop order for the three- 
point correlation C3 by the "direct steepest-descent" re- 
summation scheme. The double lines are the nonlinear two- 
point functions C and R, while the single lines are the linear 
two-point functions Cl and Rl as in Fig. [51 The black dots 
are again the three- leg vertex Kg. 

Of course, it is also possible (and more common) to 
use the standard diagrams for the computati on of the bis- 
pectrum within standard pertu rbation theory (jScoccimarrol 
Il997t IScoccimarro et al1ll998D . As noticed above, for our 
purposes the interest of the description of Fig. [2] is that 
it clarifies the link with the "direct steepest-descent" re- 
summation scheme that we also investigate in this article. 
We focus on this specific resummation scheme to be con- 
sistent with our previous stu dy for the power spectrum in 
[Valagc as fc Nishimichil ()2011[ ). In addition, as discussed in 
that article (in particular in its Appendix A), this method 
allows a fast numerical implementation. Moreover, its pre- 
dictions for the bispectrum (and higher order correlations), 
and more precisely the structure of it s diagrammatic ex- 
pansion, have already been studied in IValageasI (|2008D , so 
that no further theoretical work is required. Then, within 
this "direct steepest-descent" resummation scheme, the di- 
agrams for the bispectrum up to 1-loop order are those 



a= 6 



+ 48 



+ 48 




+ 



(f) 



+ 24 




"mixed" 



Fig. 4. Diagrams obtained at 1-loop order for the three- 
point correlation C3 within a "mixed" case. We keep the 
resummed diagram (a) of Fig. [3] but we replace diagrams 
(f), (g), (h), and (i) of Fig.|n]by their leading contributions, 
given by diagrams (f), (g), (h), and (i) of Fig. [21 



shown in Fig. [31 The difference with Fig. [51 is that we now 
have double-line two-point functions. They correspond to 
the nonlinear two-point correlation and response functions, 
as given by the same method up to the same 1-loop order. 
The single lines are the linear two-point functions, as in 

Fig.ia 

Within this resummation scheme, nonlinear two-point 
functions at "1-loop" order actually contain an infinite 
number of "bub ble" diagrams, as shown in Figs. 8 and 9 of 
IValageasI (|2008D . There, the label "1-loop" does not mean 
that we only include 1-loop diagrams, as in the standard 
perturbation theory of Fig. [51 but that the diagrammatic 
expansion is only complete up to 1-loop (i.e., although we 
include some diagrams at all orders, we miss some contribu- 
tions at 2-loop and higher orders). Then, substituting the 
expression of the nonlinear two-point functions in terms of 
the linear functions, we can check that the diagram (a) of 
Fig. [factually contains the five diagrams (a),..,(e), of Fig. [51 
as well as an infinite number of higher order diagrams. On 
the other hand, the diagrams (f),..,(i), of Fig. [31 contain 
their counterparts of Fig. [51 each diagram among (f),..,(i), 
of Fig.[51is the lowest order contribution to the correspond- 
ing diagram in Fig. [31 where nonlinear two-point functions 
are replaced by their linear counterparts. 

At the order Pi, which corresponds to the 1-loop di- 
agrams of Fig. [51 we can also consider the "mixed" case 
shown in Fig. 21 where we keep the resummed diagram (a) 
of Fig.[3l but we use for diagrams (f),..,(i), their lowest order 
terms shown in Fig. [51 Of course, these three choices, drawn 
in Figs. [51 [31 andlH agree up to order P|, and only differ 
by the number of higher order terms that are included. 

Thus, for the perturbative (3-halo) contribution (13^ 
we will consider the three alternatives shown in Figs. [51 



31 and m Contr ary to the numerical computations used in 
ValageasI (l2008l) . we do not introduce any approximation 
for the computation of these diagrams. Therefore, the per- 
turbative term i?pert(fci, ^2, ^3) contains no free parameter 
and is exact, up to 1-loop order (or to the truncation order 
of the perturbative scheme). This is an improveme nt over 
most previous studies involving the halo model (,Ma fc Frvl 
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l2000bt iFosalba et al.ll2005h , which only used the tree-order 
bispectrum for the 3-halo term, with the addition of 
bias factors (which we do not introduce in our approach) . 



4. Comparison with numerical simulations 

4.1. Numerical procedure 

We use a set of large iV-body simulations in a ACDM 
universe, consistent with the five-year observation of the 
WMAP s atellite dKoma tsu ct al. 2009), which are de- 
scribed in IValageas &: Nishimichil (|20lil) . We analyze the 
four main simulations {N = 2048'^), together with supple- 
mentary smaller simulations (N = 1024"^, 512"^ and 256'^) 
for convergence tests. 

We measure the bispectrum from a fast Fourier trans- 
formation of the density field obtained by the cloud- 
in-cells interpolation of the A^-body particles, at z = 
0.35, 1 and 3. In doing so, we use the folding scheme 
to speed up the measurements at large wavenumber 
bins without systematic error from the interpolation 
(see IValageas fc Nishimichil ()2011|) and also IColombi et al.l 
(I2009t) 1. We set the wavenumber bins for fci, ^2, and k^, 
logarithmically with 2 bins per factor 2. 

The statistical uncertainties (at l-tr level) of the mea- 
surements are estimated assuming that they main ly come 
from their Gaussian part (jScoccimarro et al1ll998[ ). 



[AB{ki,k2,k3f = 



(27r)3 2Mri 



P{k,)P{k2)P{k3), (40) 



where the factor S123 is 6 for equilateral triangles, 2 for 
isosceles triangles and 1 otherwise. In the above, V denotes 
the volume of the simulations and A^tri is the number of 
Fourier space triangles that fall into the bin of (fci, fc2, ^3)- 
We do not correct for the shot-noise contributions to the 
power spectrum and the bispectrum because the naive ex- 
pectation for the shot noise obtained by assuming a Poisson 
process does not seem very accurate, especially on large 
scales. Because we start the simulation from a regular lat- 
tice, the discreteness noise is boxed in on scales smaller than 
the inter-particle distance at the beginning. After gravi- 
tational evolution, the density fluctuations neatly follow 
the linear growth rate at fc < 0.1/iMpc~^, hence we con- 
clude that we do not have any sign of shot noise on those 
large scales. On small scales, however, we can see that our 
meas urements from sim ulations approach the Poisson noise 
(e.g., iVerde et al.ll2002h . 



Pshot(fc) 
-Bshot(fcl,fc2,fc3) 



1 



(27r)3 N 



(41) 



{2n) 



^(P(fci) + P(fc2)+P(fc3)) 



N 



(42) 



for the power spectrum and the bispectrum, respectively. 
Note that the shot noise is included in P(fci), P{k2), and 
P(fc3) of Eq. (I42p . Instead of subtracting these contribu- 
tions from the measured spectra, we assess the shot-noise 
level from Eqs. (^1]) and (02]), and then plot their relative 



importance in Figs. [T51 and [T71 below. For the reduced bis- 
pectrum, we estimate the shot noise level from 



g, 



cq, shot 



(fc) 



gcq(fc) - -Beg, shot (fc) Pcq(fc) 



3(P(fc)-Pshot) 



3F(fc)2 



(43) 



which is plotted in Fig. [18] below, where again P(k) and 
Beq(fc) are the measured power spectrum and bispectrum 
(and thus include the shot-noise contribution and other nu- 
merical errors, as compared with the theoretical values). 

We now compare our model, described in Sects. [5] and 
[3] with these numerical simulations. 

4.2. Bispectrum B{ki, k2, k^) 

4.2.1. Equilateral triangles 

We show in Fig. [S] our results for the bispectrum, from lin- 
ear to highly nonlinear scales, for equilateral configurations. 
Here we take the perturbative 3-halo contribution equal to 
the standard 1-loop result, which also corresponds to Fig. [5] 
We can see that we obtain a reasonably good agreement 
with the numerical simulations over all these scales. This 
is remarkable since our model contains no free parameters. 
Indeed, the 2-halo and 1-halo terms are fully determined 
by the halo mass function and density profiles that were 
used in IValageas &: Nishimichil ()201lD for the power spec- 
trum. This provides a significant i mprovement over t he phe- 
nomenological model presented in iPan et aD (120071 ) , based 
on a generalization of the scale transform ation introduced 
for th e two-point correlation function by iHamilton et al.l 
(|l99lf ). which breaks down on highly nonlinear scales. 

However, we can see that at high redshift (right panel 
at 2: = 3) we underestimate the bispectrum in the transi- 
tion range between linear and nonlinear scales. The same 
behavior appears for the po wer spectrum, see Fig. 9 in 
IValageas fc Nishimichil (|201ll) . This is likely to be due in 
part to higher order perturbative terms, which play a 
greater role at z = 3 than at z = 0.35, in agree ment with 
the detailed study performed in IValageasI (|201lh that com- 
pares perturbative and nonperturbative contributions. On 
the other hand, it is natural to expect the transition range 
to be the most difficult to reproduce by models of the kind 
studied in this paper. Indeed, this domain is already beyond 
perturbation theory but does not yet correspond to the in- 
ner relaxed cores of virialized halos. Therefore, it is at the 
limit of validity of the two ingredients (perturbation theory 
and halo model) used in our approach. We will discuss in 
more detail this transition range and possible improvements 
on these scales in Sects. 16.21 and iTl 

We clearly see in Fig.[5]the decay on large scales of the 1- 
halo and 2-halo contributions, in agreement with Eqs. ([55]) 
and ([33]) . As explained in Sect. [2] this is due to the new 
counterterms W(kjqM) of Eqs. ([M]) and (|3ip. which ensure 
a physically meaningful behavi or. O n the other hand, in 
agreement with Sefusatti et al.] (|2010t) , we can see that tak- 
ing into account the 1-loop perturbative contribution signif- 
icantly extends the domain of validity of the 3-halo pertur- 
bative term, as compared with the tree-level contribution. 

4.2.2. Isosceles triangles 

We show in Fig. [6] our results for isosceles triangles as a 
function of scale fc and for the fixed shapes defined by 
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Fig. 5. Bispectrum Bcq{k) = B{k,k,k) for equilateral configurations, at redshifts z = 0.35,1, and 3. The symbols are 
the results from the numerical simulations. We show the bispectrum obtained at tree order (black solid line), standard 
1-loop order (blue dot-dashed line), using the full model (red solid line), as well as the 2-halo (cyan dotted line) and 
1-halo (cyan das hed line) contributions. The magenta dot-dashed line labeled "Pan" is the phenomenological model of 
iPan et al] (l200l . 




Fig. 6. Bispectrum B{ki,k2, fca) for isosceles configurations, ki — k2 = k and k^ = k/a, at redshifts z — 0.35, 1, and 3. 
We consider the cases a = 2 (upper row) and a ~ 4 (lower row). The symbols are the same as in Fig. [5j 



ki — k2 ^ k and fca = k/a, with constant ratios a. We con- 
sider the cases a — 2 and a = 4. As expected, as a increases 
(i.e., the triangles get "squeezed"), the 2-halo contribution 
becomes more important compared with the 1-halo contri- 
bution, because the scale l/fca grows and it is more likely 
to have two distinct halos separated by such a large dis- 
tance rather than a single halo of this size. However, for 
these moderate values of a the 2-halo contribution always 
remains subdominant because the perturbative (i.e., "3- 



halo" ) contribution is still dominant when the 1-halo term 
becomes larger than the 2-halo term. For more "squeezed" 
shapes, an intermediate regime would appear, dominated 
by the 2-halo term. 

We can check that our model agrees well with the nu- 
merical simulations and is able to describe both large and 
small scales. Again, this extends the analytical predictions 
be yond the scales de scribed by the phenomenological model 
of jPan et al.1 (|2007() . However, at z = 3 (right panels) we 
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Fig. 7. Bispectrum B{ki, k2, fcs) for isosceles configurations, ki = k2, at redshifts z = 0.35, 1, and 3 (from top to bottom 
in each panel). We show our results as a function of fci = ^2 at fixed ks (upper row), and as a function of ks at fixed 
ki = k2 (lower row). The symbols are the same as in Fig. [S] 



can see the underestimation on transition scales already 
noticed in Fig. [S] 

Next, we show in Fig.[7]the evolution of the bispectrum 
with the shape of the triangle formed by {ki, k2, ka}, for 
isosceles configurations fci — k2. 

First, we consider in the upper row the dependence on 
ki = k2 = k Sit fixed fcs, so that the two equal-length sides 
run from k^/2 (flat triangle that reduces to a line) to infinity 
(squeezed limit). In agreement with the behavior observed 
in Figs. [5] and [6l the bispectrum decreases for higher values 
of fci — k2. The value of the fixed side k^ grows from the 
left to the right panel, so that we go from linear to highly 
nonlinear scales (in each case we consider the typical range 
A;3/2 < fci = ^2 < lOfcs). In the left panel, we are domi- 
nated by the 3-halo perturbative contribution, and since k^ 
remains small, for high values of fci = k2 the 2-halo term 
is larger than the 1-halo term. In the middle panel, we are 
already dominated by the 1-halo term (except at low k at 
z = 2) where the 3-halo perturbative contribution is larger) , 
and the the 2-halo term is subdominant. In the right panel 
the bispectrum is completely dominated by the 1-halo con- 
tribution. 

Second, we consider in the lower row the dependence on 
k'i at fixed value of the common sides ki = k2 = k, so that 
the third side runs from (squeezed limit) to 2ki (flat trian- 
gle that reduces to a line) . Again, the bispectrum decreases 
as the dependent wavenumber, here fca, grows. However, the 
dependence is much weaker than the one found in the upper 
row. In the left panel, we are dominated by the linear-order 



perturbative contribution, so that all curves agree with each 
other. The too high values of the numerical points at low 
fca are a numerical artifact caused by the finite boxsize and 
the numerical binning used to compute the bispectrum, be- 
cause linear theory must be increasingly accurate on larger 
scales. In the middle panel, we can see the crossover be- 
tween the 3-halo and 1-halo regimes. The 2-halo term is 
never dominant on these scales. In the right panel we are 
mostly dominated by the 1-halo term. 

In all cases we can see that our model agrees well with 
the numerical simulations. However, at 2; = 3, in the upper 
middle panel and in the lower right panel we can see a trace 
of the underestimation on the transition scales found in the 
right panels of Figs. [5] and [HI 

4.3. Reduced bispectrum 

As seen in the previous figures, the bispectrum varies 
by many orders of magnitude on the scales of interest. 
Therefore, it is customary to introduce the "reduced bis- 
pectrum" defined as 



Q{ki,k2,ks) 



B{ki, k2, fcs 



P(fci)P(fc2)+2cyc. 



(44) 



Indeed, this ratio goes to a constant on large scales, as 
seen from the tree-order express ion (1391). while it shows a 
moderate growth on small scales (jBernardeau et al1l2Q02aD . 

We compare in Fig.[8]our results for the reduced equilat- 
eral bispectrum, Qeq(fc) — Q{k, k, k), with numerical simu- 
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Fig. 8. Reduced bispectrum, Qcq{k) = i?cq(fc)/[3P(fc)^], for equilateral configurations, at redshifts z ~ 0.35,1, and 3. 
The symbols are the same as in Fig. [5] The vertical arrow in the upper right part shows the wavenumber beyond which 
the simulation shot noise is greater than 10%. 



lations. Here, in the denominator of Eq. p^ we use for the 
theoretical predictions the power spec trum ob tained with 
our model described in Valageas & Nis himichil (|2011) (that 
is, using the same halo model as in this paper and the di- 
rect steepest-descent resummation for the perturbative 2- 
halo term), while for the numerical curves we use the power 
spectrum measured in the simulations. 

The vertical arrow in the upper right part of each panel 
of Fig. [S] shows the wavenumber where the shot-noise level 
of the numerical simulations becomes larger than 10%. We 
can see that this also roughly corresponds to the scale where 
the measured reduced bispectrum is maximum. Therefore, 
the downturn and the decrease of the reduced bispectrum 
at higher k found in the simulations is only a numerical 
artifact caused by the limited resolution, and these points 
must be discarded. In particular, at lower redshift where 
the simulations are able to probe deeper into the nonlinear 
regime, we clearly see in the left panel at z = 0.35 the 
fast growth of Qeq{k) in the highly nonlinear regime, after 
a small plateau on the transition scales. As shown by the 
two left panels, this high-fc fast increase of Qcq{k) is well 
reproduced by our model. 

We can see that we obtain a reasonably good agreement 
with the simulations on both large and small scales. This 
is consistent with our previous results, shown in Fig. [S] and 
in|Valageas & Nishimichi (2011), where we found that both 
the bispectrum and the power spectrum are well reproduced 
on quasilinear and highly nonlinear scales. As expected, we 
recover the large-scale plateau associated with the tree-level 
limit pop , and the early rise owing to higher order pertur- 
bative contributi ons. At very high k, in agreement with 
previous studies (iMa k Frvl l2000bl: ICoorav fc Shetbl l2002t 
iFosalba et al.ll2005[ ). the halo model leads to a continuous 
growth of the reduced bispectrum, which is a marked dif- 
ference with the prediction of the stable clustering ansatz 
lPeeblesl (,1980'). This agrees with our numerical simulations 
u p to the scales th at are well resolved. As shown in Sect. 2 
of lValageasI (|1999D . the positivity of the matter density ac- 
tually implies that the reduced bispectrum, and more gen- 
erally the real-space coefficients Sp — {p^)c/{pji)c~^j reach 
a constant or keep growing on small scales, as the den- 
sity variance (/9fl.)c increases. The limiting behavior of the 
constant reduced bispectrum and the constant coefficients 
Sp is achieved for the stable clustering ansatz, and more 



generally for multifractal models such that (p^)c are gov- 
erned by a single fractal exponent a for the values of p that 
are considered. The density field described by a halo model 
clearly violates these condition^ because it does not display 
this scale invariance, with a characteristic nonlinear mass 
associated with the falloff of the halo mass function and 
reasonably smooth profiles that depend on the mass scale 
(through their concentration parameter) . This implies that 
Qcq{k) has to grow on small scales, as checked in Fig. |S1 

On the transition scales, we underestimate both the bis- 
pectrum and the power spectrum, and it appears that the 
latter effect is dominant, which leads to the "bump" seen 
in Fig. [51 As seen from the numerical results, this feature 
is not physical and only caused by the shortcomings of 
the theoretical model on this ra nge. Such a feature was 
also noticed in previous studies (IScoccimarro et al.l 120011 : 
iTakada fc Jainll20 03': 'Fosalba et al "2005V which found that 
by using a larger halo radius, taking into account exclusion 
constraints in the 2-halo and 3-halo terms, or introducing 
a high-mass cutoff in the halo mass function, one might re- 
move this "bump" and reach a better agreement w i th nu- 
merical simulations. On the other hand. lGuo fc Jin3 ()2009() 
find that no self-consistent halo model is able to reproduce 
well both the matter power spectrum and bispectrum on 
these transition scales. We will return to this problem and 
these modifications in Sect. 16.21 and will discuss in Sect. [7] 
an alternative procedure to improve the model on these 
scales, using this artificial "bump" as a signature of insuf- 
ficient accuracy and introducing simple interpolations that 
resolve most of this problem, while keeping large and small 
scales unchanged. 

5. Comparison with previous halo models 

We show in Fig. [S] the impact of two features of our model 
that differ from previous implementations of the halo model 

^ The halo model can be made to recover the stable-clustering 
ansatz predictions if t he mass function scales at low m ass as 
n{M)dM oc dM/M^ (|Valageasl [19991 : iMa fc Frvll2000bD . This 
unrealistic formal limit, where the apparent amount of matter 
per unit volume is infinite, corresponds to a multiple counting 
of "halos", which contain an infinite hierarchy of substructures 
that are also counted in the mass function, in agreement with a 
fractal model. 
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Fig. 9. Ratio B{ki, k2, k3)/[Ps{ki)Ps{k2) + 2cyc.], where Ps{k) is the fit to the power spectrum froml Smith eTall (pOOl . 
The red sohd line is our model, as in Figs. [5] to [51 the black dotted line is the approximation Ban = Btree, the green 
dot-dashed line corresponds to no counterterms in the 1-halo and 2-halo contributions {W — 0), and the blue dashed 
line makes both approximations. Upper row: equilateral configurations, at redshifts z = 0.35, 1, and 3, as in Fig.[SJ Lower 
row: various isosceles configurations at redshift z = 0.35, as in left panels of Figs. Eland [T] 



(IScoccimarro et al.|[200lt iFosalba etall 120051: iGuo fc Jind 
I2009D . i) the use of 1-loop perturbation theory instead of 
the tree-level contribution ((39)) for the 3-halo term, and ii) 
the counterterms W in Eqs. ([M)) and (PT|) in the 1-halo and 
2-halo terms. 

To distinguish more clearly between various models, 
that is, to reduce the size of the vertical axis, we plot 
in Fig. ini the effective reduced bispectrum defined as in 
Eg. (1311) ■ but where we use the power spectrum Ps(k) from 
ISmith et al.l (|2003l ) in the denominator, for the models as 
well as for the numerical simulations. Thus, by dividing the 
bispectra by the same denominator we avoid introducing 
additional errors through the denominator and can make 
a clear comparison between the various bispectra and the 
simulations. In other words. Fig. ^ is not meant to show 
the actual reduced bispectrum (|44| , which was displayed in 
Fig-Hlfor equilateral triangles, but only to show the bispec- 
trum i? on a more convenient scale. We focus on quasilin- 
ear scales where the different models can be distinguished, 
since on very large or small scales they converge towards 
the tree-level contribution (|M1) or the 1-halo contribution 
with a negligible counterterm. 

First, we can see in Fig. |9l that discarding the con- 
tribution of 1-loop diagrams (the diagrams (b) to (i) in 
the representation of Fig. ^ leads to a significant un- 
derestimate of the bispectrum on weakly nonlinear scales, 
k ~ 0.2/iMpc^i at z < 3, in agreement with previous works 



based on perturbation th eory alone (|Scoccimarro et all 
119981: ISefusatti et"aI1l2010[ ). In particular, at z = 0.35 the 
1-loop contribution is necessary to obtain a good match to 
the simulations and appears to be sufficient to make the 
bridge to the smaller scales where the 1-halo term is dom- 
inant. At z = 3, the 1-loop contribution also extends up 
to fc ~ 0.4/iMpc~^ the good agreement with simulations, 
while using the tree-level contribution alone yields signifi- 
cant discrepancies as soon as fc > 0.2/iMpc~^. However, it 
is no longer sufficient to make the bridge with the highly 
nonlinear scales dominated by the 1-halo term, because we 
now underestimate the bispectrum on the transition scales, 
k ~ l/iMpc"^. This behavior was already noticed in Figs. [SI 
andini 

Second, setting the counterterms to zero in the 2-halo 
and 1-halo contributions (while keeping the 1-loop contri- 
bution in the 3-halo term) obviously increases the bispec- 
trum and worsens the agreement with simulations. This is 
most clearly seen on the quasilinear scales, which are al- 
ready well described by standard perturbation theory so 
that the extra power associated with the unphysical con- 
stant asymptotes at low k of the uncorrected 2-halo and 
1-halo contributions spoil the good agreement with simula- 
tions. 

Making simultaneously both approximations, discard- 
ing 1-loop diagrams and halo counterterms, as in usual im- 
plementations of the halo model, also yields significantly 
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worse results. Although both effects partly compensate, ne- 
glecting 1-loop contributions dominates so that the result- 
ing bispectrum is significantly too low on weakly nonlinear 
scales. Therefore, in addition to a more satisfactory phys- 
ical behavior (systematic agreement on large scales, up to 
higher order, and decay of 2-halo and 1-halo contributions), 
our approach gives a better accuracy on weakly nonlinear 
scales. This shows the importance of including both 1-loop 
contributions and the 2-halo and 1-halo counterterms. 

6. Dependence on various ingredients of the model 

We now investigate the impact of various ingredients of the 
model on the predictions obtained for the bispectrum. 

6.1. Dependence on the perturbative scheme 

We investigate in this section the dependence of our results 
on the perturbative scheme used for the 3-halo contribu- 
tion. As in Fig. [HI the bispectra obtained from the analyt- 
ical models and the numeric al simulations are d ivided by 
the same power spectra, from lSmith et al] (|2003[ ). to make 
a clear comparison. 

We compare in Fig.[TU]the results obtained using for the 
perturbative 3-halo contribution either the standard 1-loop 
result of Fig. [5J the full steepest-descent resummation of 
Fig. [21 or the mixed case of Fig. HI For each choice we show 
both the 3-halo term itself (blue curves that decrease at 
high k) and the full prediction, where we sum the 3-halo, 
2-halo, and 1-halo contributions. (The 2-halo and 1-halo 
contributions are the same for all three choices.) 

Let us first consider the upper row, which corresponds 
to equilateral triangles. (At 2 = 3 we can see the underesti- 
mate of the bispectrum in the transition range that we had 
already noticed in Fig. [S]) The standard and mixed cases 
give very close results. At z = 3, the mixed case, giving a 
slightly greater 3-halo contribution, yields a slightly better 
agreement with numerical simulations around l/iMpc~^. 
However, this is not conclusive, especially since the 1-halo 
term is already non-negligible on these scales and it is not 
sufficient to significantly improve the match to simulations. 
The complete steepest-descent resummation of Fig.[2]yields 
a 3-halo contribution that decays significantly faster at high 
fe. At z = 0.35 and z — 1, two numerical points in the range 
0.1 < k < 0.2/iMpc~"'^ suggest that this may be a true im- 
provement on these quasilinear scales. This improvement 
might be important for the analysis of baryon acoustic os- 
cillations; it has been shown that resummation schemes 
give bette r predictions for the power s pectrum on this 
scale (e.g., 'Crocce & Scoccimarro' (2008[): ]Nis"himichi et ahl 
(|2009l) : [Valageas & Nishimichi (2011), and we could ex- 
pect the same thing for the bispectrum. However, beyond 
0.2/iMpc~^ the faster decay significantly worsens the agree- 
ment with the simulations. 

The upper row of Fig.[TU]does not allow us to clearly dis- 
criminate between the various schemes because they mainly 
differ in the transition regime, where the 1-halo contribu- 
tion is already important and all models tend to underes- 
timate the bispectrum. However, looking now at the lower 
row, associated with isosceles triangles, and especially at 
the two lower right panels that show the evolution of the 
bispectrum with the shape of the triangle, we can see that 
the standard 1-loop perturbation theory yields a signif- 
icantly better agreement with simulations than the two 



other perturbative expansions. This is also the simplest 
scheme for practical purposes. 

Thus, from the analysis of both equilateral and isosce- 
les triangles, we can conclude that using the standard 1- 
loop perturbation theory for the perturbative 3-halo term 
is the most efficient and accurate scheme among the three 
approaches studied in this paper. This is quite different 
from t he power spectrum, studied in I Valageas fc Nishimichil 
(|201l[ ). where the 1-loop steepest-descent resummation is 
clearly better than standard perturbation theory. There, it 
is more accurate on weakly nonlinear scales but it is also 
the only choice (with other resummation schemes) that al- 
lows a combined model for both large and small scales. 
Indeed, the standard 1-loop perturbation theory yields a 
2-halo contribution that keeps growing on small scales (for 
the logarithmic power of Ea. (j45p below) and that worsens 
the agreement with simulations (even though it is smaller 
than the 1-halo term). 

We do not have this small-scale problem for the bis- 
pectrum, as seen in Figs. [SJ El and [Hi Indeed, the 1-loop 
standard perturbation theory gives a contribution to the 
bispectrum that quickly decreases at high k and does not 
spoil the good agreement with simulations shown by the 1- 
halo term in this regime. This is best seen in Fig. |8l which 
shows that the 1-loop standard perturbation theory contri- 
bution to the reduced equilateral bispectrum Qcq{k) itself 
also decreases at high fc, in spite of the denominator 3P(fc)^. 
Therefore, contrary to the power spectrum, it is now possi- 
ble to use the 1-loop standard perturbation theory to build 
a combined model that covers all scales. Moreover, as shown 
by Fig. [TUl it happens that this is also more accurate than 
the two other perturbative schemes investigated in this pa- 
per, which we presented in Figs. [3]and|4l 

It is possible that other perturbative schemes that we 
did not study here provide more accurate results. On the 
other hand, higher orders of standard perturbation theory 
may yield contributions that are increasingly large on small 
scales, so that beyond a certain order they lead to a 3-halo 
term that is unphysically large and spoils the good agree- 
ment with simulations on small scales. Then, one would 
need to use other perturbative approaches, such as the re- 
summation schemes investigated here, or to add some ad- 
hoc cutoff on small scales. 

For completeness, we mention a few differ- 

ent app roaches. A first method, introduced in 

ICrocce fc Scoccir narrol ()2006bl lal) from a diagrammatic 
approach, reorganizes the standard perturbation theory 
and performs partial resummations by introducing both 
correlation functions and propagators (or response func- 
tions), in a fashion somewhat similar to the approach 
described in Sect. 13.21 However, a key step in this method 
is the matching between the low-order low-fc behavior 
and the high-A: asymptote of the propagator. This en- 
sures a good behavior of this quantity in the nonlinear 
regime, whi ch has been checke d again s t nu merical simu- 
lations (Crocce fc Scoccimarrol l2006al : iBernardeau et al.l 
200^, and expresses a well- understood "sweeping ef- 
fect" associated with the random tran sport of density 
structures by large-scale velocity flows ()Valageasl l2007bl : 
Bcrnar deau fc Val ageas '2008'. '2010\ Another advantage of 
this approach is that its e xtensions to high-order quantities 



(IBernardeau et aL 20081) . to non-Gaussian initial condi- 
tions (^Bernardeau et al."2 010t) . and to higher perturbative 
orders (Anselmi et al... 20101) , have already been studied. 




Fig. 10. Ratio B{ki, k2, k^) / [Ps{ki) Ps{k2) + 2cyc.], as in Fig. [S] We show the resuhs obtained by using for the 3-halo 
term each of the three 1-loop perturbative expansions given in Figs. [5J [31 and S) In each case, we show both the 3-halo 
contribution alone (which decays at high fc, upper line in front of the labels) and the full result (i.e., the sum of the 
3-halo, 2-halo, and 1-halo terms, lower line in front of the labels). We also plot the tree-level perturbative result (black 
solid line). Upper row: equilateral configurations, at redshifts z = 0.35, 1, and 3, as in Fig. [SJ Lower row: various isosceles 
configurations at redshift z — 0.35, as in left panels of Figs. [S] and [T] 



A second approach, developed in IValageai ( 2007allbl) as 



a "larg e-A^ 2PI" exp ansion, and in iTaruva fc Hiramatsul 
(pOOSi) : iTaruva et all (pOOl as a "closure approximation", 
is quite similar to the one described in Sect. 13.21 Although 
in principle it may give more accurate results, its full im- 
plementation is more complex because self-energy terms 
depend on the nonlinear quantities that are looked for. 
This leads to coupled nonlinear equations, which depend 
on scale s and times, an d to heavier numerical compu- 
tations (jValageasI l2007aD . Because for practical purposes 
one requires fast algorithms this presents a significant dis- 
ad vantage, unless one introduces further simplifications 
(jTaruvaet al.ll2009t) . 



A third method, introduced in lPietro"nil ()2008[ ). directly 
obtains the hierarchy of equations obeyed by the many- 
body correlation functions from the equations of motion, in 
a fashion similar to the standard BBGKY hierarc hy. Then, 
trunca tion at a given order (at the trispectrum in iPietronil 
()2008[) ) defines an approximation for lower-order correla- 
tions. This method ha s already been used to stu dy the effect 
of massive neutrinos (Lesgourgues et al. 2009) and of pri- 
mordial non-Gaussianities (Bartolo et al. 2010 ). An advan- 
tage of this approach is that it does not involve propagators, 
but only the usual many-body correlations also encountered 
in standard perturbation theory, and it is written in terms 



of single-time quantities. This is a great simplification that 
should allow efficient numerical computations. 

A fourth approach, developed in iMatsubaral (|2008bllah . 
is based on a Lagrangian framework. This would be 
well-suited to the approach described in this paper and 
presents the advantage that it provides direct extensions 
to redshift-space stati st ics. U nfortunately, as noticed in 
IValageas fc Nishimichil ()2011[ ). in its present form this 
Lagrangian perturbation theory does not fare as well as its 
Eulerian counterparts when tested against numerical sim- 
ulations. However, this approach remains interesting, espe- 
cially in view of applications to redshift space. 

In any case, it would not be too difficult to combine any 
of these methods with halo models by using the framework 
described in this paper, which is more general than the use 
of the standard perturbation theory or the two resumma- 
tion schemes described in Sect. 13.21 



6.2. Dependence on halo properties 

We now investigate the dependence of our results on the 
details of the halo model. As we recalled in Sect. 14.31 
the "bump" shown by the reduced bispectrum in Fig. [5] 
was also found in previous works, which noticed that 
this feature ma y be cured to some e xtent by tu ning halo 
parameters (,Scoccimarro et al.l 120011 : iTakada fc'jain 20031 
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Fig. 11. Reduced bispectrum (upper panel, as in Fig. [S]) 

and the bispectrum scaled by 3Ps(fc)^ (lower panel, as in the 
upper row of Fig. ITU)) obtained for modified halo properties 
on equilateral triangles at z = 0.35. We plot our fiducial 
model (red solid line, as in previous figures), the impact 
of a high mass cutoff (M < l{)^^h~^MQ, blue dashed line; 
M < 5 X lO^'*/i~^M0, green dot-dashed line), and the effect 
of defining halos by the lower nonlinear density threshold 
(5 = 50 (black dotted line). 



iFosalba et al]l2005D . Therefore, we show in Fig. [TT] the re- 
duced bispectrum (upper panel, as in Fig. and the scaled 
bispectrum (lower panel, as in Fig. llOp for several modifi- 
cations of halo parameters at redshift z = 0.35. 

FollowinglScoccimarro et all (|200lf ): IWang eraII(|2004D 



IFosalba et al 



(120051) . we investigate the impact of trun- 
cating the halo mass function below M < 10^^ or M < 
5 X 10^'^h~^ Mq. In agreement with these studies, removing 
large halos decreases both the matter power spectrum and 
bispectrum on mildly nonlinear scales (since smaller scales 
are associated with small halos), and the net effect on the 
reduced bispectrum is also a small decrease of the artificial 
"bump". However, the upper panel shows that even trun- 
cating at M < 5 X 10^^/i~^Mq is not sufficient to erase 
the "bump" and to provide a good agreement with simula- 
tions for the reduced bispectrum. Two differences with pre- 
vious studies are that the 3-halo term includes 1-loop per- 
turbative contributions, which are still important on these 
scales as shown in Figs. 151 and ITUl and that the 1-halo and 
2-halo contributions contain the counterterms W, so that 
the relative importance of the 2-halo and 1-halo contribu- 
tions is somewhat smaller than in previous models on the 
weakly nonlinear scales associated with the early rise of the 



Fig. 12. Reduced bispectrum (upper panel) and the bispec- 
trum scaled by 3Ps{k)^ (lower panel) obtained for modified 
halo properties on equilateral triangles, as in Fig.[Tl]but at 
z = 3. 



"bump" . On the other hand, because both the power spec- 
trum and the bispectrum simultaneously decrease when we 
add such a high-mass cutoff and these effects partly com- 
pensate in the reduced bispectrum, before we obtain a sig- 
nificant improvement for the latter quantity both the power 
spectrum and the bispectrum have already been decreased 
by a large amount that disagrees with numerical results, as 
shown in the lower panel for the bispectrum. Moreover, in 
our large simulations we find halos up to 3 x 10^^h~^ Mq, 
so that the cutoff at M < 5 x 10^'^h~^MQ is already too low 
to be fully justified by the finite box size of the simulations. 

Then, following IFosalba et al.l ()2005D . we investigate the 
impact of larger halo radii. More precisely, we consider the 
impact of defining halos by a nonlinear density threshold 
S — 50 instead of 6 = 200, following the p r ocedu re de- 
scribed in Sect. 6.1 of I Valageas fc Nishimichil (|201lD . This 
also involves a slightly different halo mass - concentration 
parameter relation, so as to keep a satisfactory match to nu- 
merical results for the power spectrum. As seen in Fig. Illl 
this modification does not bring a significant improvement 
either. 

We find similar results at higher redshifts, as seen in 
Fig- [m at z = 3 (where we consider the smaller mass thresh- 
olds Af < 3 X 10" and M < lO^^h-^MQ). 

Th erefore, we reach the same conclusion as lGuo &: Jinei 
()2009D . that tuning these halo parameters cannot simul- 
taneously provide accurate results for the power spectrum 
and bispectrum, especially compared with large simulations 
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where massive halos are present, which removes the justifi- 
cation for introducing severe high-mass cutoffs. 

The discrepancies found on these transition scales are 
not so surprising because one expects the transition range 
to be the most difficult to describe in a systematic fashion, 
since both the perturbative expansion (that neglects shell 
crossing) and the halo model (which assumes spherical and 
relaxed objects in our implementation) may break down in 
this intermediate regime. Another possibility is that we are 
missing important high-order perturbative terms that have 
not been included in the perturbative expansions used here, 
which are only complete up to one-loop order. Indeed, for 
the power spectrum perturbative terms up to ord er ^ 9 (at 
z = ) or ^ 66 (at z = 3) are likely to be relevant (jValageasI 
l201lt) . On the halo-model side, because in reality the den- 
sity field on these transition scales shows a crossover from 
relaxed inner halo shells to outer infalling shells and fila- 
ments that cannot any longer be described in terms of well- 
defined halos, one cannot expect a systematic convergence 
to the numerical results by adding such simple modifica- 
tions to halo properties. Then, one is probably limited to 
some extent by the intrinsic limitations of the halo model 
itself. 

7. Improving the predictions for P{k) and Bcq{k) 
using the reduced bispectrum Qeq{k) 

The artificial "bump" found in Fig. |S] for the reduced equi- 
lateral bispectrum Qcq{k), which is mostly caused by the 
underestimation of the power P{k), suggests a simple trick 
to improve the predictions for P{k) and Bcq{k). The idea is 
to use the unphysical "bump" shown by the predicted re- 
duced bispectrum Qeqik) to automatically detect the range 
where the model is not sufficiently accurate. The 
procedure that we investigate in this paper is to draw in the 
(log fc, Qcq) plane the lower tangent line to the predicted 
curve that was plotted in Fig. [51 We show this construction 
in Fig. [131 where we again plot the prediction of our model, 
described in the previous sections and labeled "direct" in 
this figure, as well as the lower tangent on the region where 
the "bump" appears. The two contact points between the 
model curve and its tangent line define the two wavenum- 
bers, k- and between which the artificial "bump" arises 
and the model needs to be improved. 

First, we modify the density power spectrum as shown 
in Fig. 1141 where we plot the power per logarithmic interval 
of k, defined as 

A'^{k) =ATrk^P{k). (45) 

The modified power "ta ng." is obtained from the "dir ect" 
prediction described in IValageas Sz Nishimichi! ()2011h by 
drawing in the (logfc,logA^) plane the upper tangent line 
on the interval that runs through the left point 

(logfc_,log A^(fc_)). For the cases shown in Fig. [Til the 
right contact point has an abscissa fc^ < fc+ (by construc- 
tion we have fc_ < A;^ < fc+), and the power spectrum is 
only modified on the interval [fc_, fc^]. As seen in Fig. [TH 
in this fashion we correct most of the artificial "dip" shown 
by our model without modifying the large-scale and small- 
scale regimes where the "direct" predictions are satisfac- 
torj0. 

^ Although we show in Fig. [14] this geometrical construction 
for the logarithmic power A^(fc), applying the same construction 
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Fig. 14. Power per logarithmic interval of fc, as defined in 
Ea.(|45|). at redshifts z = 0.35, 1, and 3 (from top to bot- 
tom). We s how the "direct" model (red solid line), stud- 
ied in Vala geas fc Nishimichil (|2011l) . and the geometrical 
modification "tang." (green dot-dashed line), given by the 
upper tangent that runs through the point of abscissa fc_ . 
The small dashed vertical line shows the location fc^ of the 
contact point of this tangent line with the "direct" curve. 
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Fig. 15. Equilateral bispectrum at redshifts z = 0.35, 1, 
and 3 (from top to bottom). We show the "direct" model 
(red solid line), already shown in Fig.[5l and the geometrical 
modification "tang." (green dot-dashed line) given by the 
upper tangent that runs through the point of abscissa fc_ . 
The small dashed vertical line shows the location fc" of the 
contact point of this tangent line with the "direct" curve. 



Second, for equilateral triangles we modify the bispec- 
trum in the same manner. Thus, as shown in Fig. 1151 we 
again draw in the (log k, log -Bcq) plane the upper tangent 
line, on the interval [/c_,fc_|_], that runs through the left 
point (log/c_,logi3oq(fc_)). This yields another right con- 
tact point at fc" , which in the cases shown in Fig. [151 is 
located to the left of k^ as there is again a knee in the 



to the power P{k), that is, in the (logfc,logP) plane, yields the 
same results (since the abscissa k'j^ of the contact point with the 
upper tangent line is the same). 
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Fig. 13. Reduced bispectrum, Qcq{k) = Bcq{k) /[3P{k)'^], for equilateral configurations, at redshifts z — 0.35, 1, and 3. 
The points are the results from the numerical simulations and the solid line is our "direct" model, as in Fig. [51 The black 
dashed line is the lower tangent to the "direct" curve Qcq{k), on the transition region. It defines the two contact points 
of abscissa k- and k+, shown by the two small vertical lines, between which we introduce the two modifications "tang." 
(green dot-dashed line) and "flat" (blue dashed line), explained in the main text. 



shape of the bispectrum. This ensures that we follow the 
"direct" prediction beyond the knee, where it shows a good 
match to the numerical simulations, while correcting most 
of the artificial "dip" shown by the model. 

As explained above, the interest of this procedure, based 
on the reduced equilateral bispectrum Qcq{k), is to auto- 
matically define the wavenumbers fc_ and fc+. One could 
imagine defining these boundaries in a simpler manner from 
the power spectrum itself, for instance by A^(A;±) = Aj_, 
with fixed values Aj_ that would mark the transition regime 
(such as = 1 and A^ = 30). However, owing to the 
curvature of the linear CDM power spectrum, the inter- 
val [/c_,fc-(_] where the model departs from the numerical 
simulations is not given by these redshift-independent con- 
ditions, as can be seen in Fig. 1141 In particular, at higher 
redshift (i.e., lower value of the effective slope n of the lin- 
ear power spectrum on the relevant scales) the threshold 
A?_ is seen to decrease. This is partly captured by our pro- 
cedure, described in Fig. [121 which is sensitive to the shape 
of the initial conditions, through the shape of the predicted 
reduced bispectrum Qoq(fc)- 

Finally, going back to Fig. [T21 from the improved 
model "tang." constructed for the power spectrum and 
the bispectrum, shown in Figs. [T31 and [121 we can build 
a new prediction "tang." for the reduced bispectrum from 
Eq. pi)) . For equilateral triangles this reads as Ql^^ {k) — 
-Seq''^ (fc)/[3F*^"s (fc)2]. The result, plotted in Fig. [H 
shows a significant improvement over the "direct" predic- 
tion, which was already plotted in Fig. [51 However, we can 
see that although it is much smaller, the artificial "bump" 
has not been fully removed by the modifications to the 
power spectrum and the bispectrum. If one is interested 
in the reduced bispectrum Qcq{k), one can introduce a 
last improvement by replacing this small "bump" by a flat 
plateau. This "flat" model is obtained from the "tang." 
curve by running down the Qeq"^' (^) curve over the interval 
[k- , k+] , starting from the right boundary k+ , and impos- 
ing a monotonic decrease as k decreases. Thus, Q^q* and 
Qoq"*' s-re identical over most of (and on all larger 

and smaller scales), except under the remaining "bump" 



of Qcq"^', where Qoq' is constant and equal to the local 
minimum of Q*™^' to the right of this "bump" . 

The geometrical improvements described in Figs. [T51[Hl 
and[T5lare not as elegant as one would wish for. Indeed, by 
combining pertur bative expansions with halo m odels, as in 
this article and in I Valageas fc Nishimichi! (|2011[ ). one would 
hope to build a model that shows a good match to numeri- 
cal simulations on all scales, without further modiflcations. 
Unfortunately, as noticed above, at the current stage there 
remain some discrepancies on the transition scales. As dis- 
cussed in Sect. 16.21 this is not surprising because transition 
scales may be at the limit of validity of both perturbation 
theory and halo models. Indeed, shell crossing (which is 
beyond the reach of the perturbative approaches studied 
here, based on the fluid approximation) is already impor- 
tant, and these scales do not correspond to relaxed halo 
inner shells, but rather to outer infalling clumps and to fil- 
aments (which cannot be described as isolated objects with 
a well-defined profile). In particular, we have seen in Fig. 1111 
that tuning halo parameters does not easily provide a signif- 
icantly better agreement with simulations. Moreover, this 
would require some ad-hoc modifications and new free pa- 
rameters, which make the model less predictive. Then, the 
geometrical modifications introduced above can be seen as 
a simple procedure to obtain a good matching between the 
weakly and highly nonlinear regimes, and are no less sat- 
isfactory than more "algebraic" matchings. They offer the 
advantage to bypass all these complicating matters and the 
free parameters they would involve, but they share with ap- 
proaches based on modifications to the halo parameters the 
lack of a systematic method toward increasingly high accu- 
racy. 

In addition to the automatic definition of the interval 
k^], associated with the transition range, an important 
property of this procedure is that the power spectrum and 
the equilateral bispectrum are not modified outside of this 
interval. Therefore, we keep the good match obtained on 
larger and smaller scales. In particular, large scales are still 
obtained by systematic perturbative expansions and small 
scales by phenomenological halo models, so that the final 
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result is still a combination of these two approaches and 
keeps their distinct benefits. 

We describe in App. |B] the impact on the real-space 
two-point correlation function of this modification to the 
power spectrum. In particular, we check that we keep an 
accuracy on the order of 1% on weakly nonlinear scales, 
X > lOft-^^Mpc, while reaching an accuracy on the order of 
10% on nonlinear scales. 

8. Typical accuracy of combined models 

We now consider the accuracy of the combined 
model described in the previous sections and in 
IValageas fc Nishimichil ()2011[ ). We first plot in Fig. [H] the 
relative difference between our model and the numerical 
simulations for the power spectrum, 

AP,,, |Pmodol(^) - -PN-body(fc)| 
P -PN-body(«) 

T he "direct" curve correspond s to the model described 
in IValageas &: Nishimichil ()2011l ) , without the geometrical 
modifications introduced in the previous Sect. [71 and was a l- 
ready displayed in Fig. 22 of lValaeeas fc Nishimichil (j201ll) . 
We can see that the modification P*™s(/j)^ shown in Fig. [HI 
in terms of the logarithmic power A^(fc), provides a signif- 
icant improvement on the transition scales, especially for 
z < 1. Thus, we obtain an accuracy on the order of 10% 
on the transition scales, and 1% on weakly nonlinear scales. 
On very large and very small scales the curves in Fig. [TBI are 
dominated by the error bars of the numerical simulations. 

Next, we show in Fig.jTTjthe relative difference between 
our model and the numerical simulations for the equilateral 
bispectrum, 

ABeq^^^ _ |-Bmodel(fc, fc, k) - i?N-body(fc, k, fc) ^^^^ 
Bcq i3N-body(fc, fc, fc) 

We can see that we reach a typical accuracy of 10% on most 
scales. At low k we are dominated by the error introduced 
by the finite size of the simulation box, as shown by the 
rise of the statistical error, and the theoretical predictions 
are actually more accurate than appears in the figure (they 
actually become increasingly good on larger scales where 1- 
loop perturbation theory is increasingly accurate) . At high 
k we are dominated by the error from the finite resolution 
of the simulations. As f or the power spectrum, studied in 
IValageas fc Nishimichil (j201lt) and shown in Fig. [161 the 
accuracy of the "direct" model is worst on transition scales, 
especially at z = 3. This agrees with the behavior found in 
Figs. [5] and [H The modified bispectrum Bl^^^ {k), shown in 
Fig. I15[ provides a modest improvement on the transition 
scales, except at z = 0.35 where in our case it happens 
to give a slightly larger error that remains on the order 
of 10%, which is typical in this range, so that the change 
caused by the geometrical modification to the bispectrum 
is not meaningful here. 

Although this 10% accuracy (and better on very large 
scales) is already a satisfactory result, it is significantly be- 
low what can be achieved for the power spectrum, where 
an accuracy on the order of 1% is reached in the weakly 
nonlinear regime, as seen in Fig. 1161 From the compari- 
son between different resummation schemes discussed in 
Sect. 16. 1[ it is not clear whether going to higher orders of 



perturbation theory would provide a significant improve- 
ment. It may happen that on these scales, especially in the 
transition regime, the bridge between the perturbative 3- 
halo term and the non-perturbative 2-halo and 1-halo terms 
is the main source of error because of the intrinsic limita- 
tions of a description in terms of relaxed spherical halos. 

On the other hand, the error bars of the simulations are 
much larger for the bispectrum than for the power spec- 
trum, and at fc < 0.1hMpc~^ the level of 10% seen in 
Fig. [T7| is mostly set by the simulations. The comparison 
with Fig. [T6| suggests that on these large scales the accu- 
racy of the theoretical model is actually much better, on 
the order of 1%, because it is determined by the systematic 
perturbation theory. Thus, our model appears to be com- 
petitive with numerical simulations because its fares as well 
or better on both large and small scales, but worse on the 
intermediate mildly nonlinear scales. Of course, for practi- 
cal purposes, the main advantage of the analytical model is 
that it allows much faster computations, as well as a greater 
fiexibility. 

Finally, we plot in Fig.[TH|the relative difference between 
our model and the numerical simulations for the equilateral 
reduced bispectrum, (A(5cq/Qcq)(fc), defined in a fashion 
similar to Egs. pB)) and (jTT)) . As in Fig. [SI the vertical arrow 
in the upper right part shows the wavenumber where the 
shot noise of the simulations becomes large. 

As for the bispectrum, we obtain a typical accuracy on 
the order of 10% (and better on larger scales), except on the 
transition scales where the prediction given by the "direct" 
model can be larger than the numerical results by up to a 
factor two. This corresponds to the artificial "bump" found 
in Fig. [5] The modified predictions Ql^'^' and Q^^* allow us 
to recover an accuracy on the order of 10% in this transition 
regime. On smaller scales the error bars of the numerical 
simulations are too large to obtain a precise estimate of the 
accuracy of these models because the rise seen in Fig. [TBlis 
due to the finite resolution. 



9. Conclusion 

Extending a previous work dedicated to the matter power 
spectrum, we have explained how to combine perturbation 
theories with halo models to build unified models that can 
describe all scales, from large linear scales to small highly 
nonlinear scales. Starting again from a Lagrangian point 
of view, instead of the usual Eulerian point of view, we 
have shown how to recover the decomposition into 3-halo, 
2-halo, and 1-halo contributions, which we relate to per- 
turbative and non-perturbative terms. This explains how 
one can build a model that agrees with perturbation the- 
ory at all orders because the 1-halo and 2-halo terms are 
non-perturbative corrections that vanish at all orders over 
Pl- Moreover, we explained how new counterterms appear 
in the 1-halo and 2-halo contributions. This ensures that 
these contributions vanish at low k, as required by physical 
constraints on the power generated by small-scale redis- 
tributions of matter. This improves previous models that 
displayed an unphysical constant asymptote at low fc, and 
allows us to reach a higher accuracy. In addition to standard 
perturbation theory, we described two alternative pertur- 
bative schemes, also complete up to 1-loop order, which can 
be used for the perturbative 3-halo contribution. They con- 
tain infinite partial resummations of higher order diagrams. 
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Fig. 16. Accuracy of our models and our numerical simulations at redshifts z — 0.35, 1, and 3 for the power spec- 
tr um. The re d solid line "direc t" is the relative difference P5|) between the simulations and our model as described 
in IValageas &: Nis himichii (|201lh . The green dashed line "tang." corresponds to the geometrical modification shown in 
Fig. [Ml The black dashed line is the error level of the simulations, including both the relative statistical error (which 
grows at low fc) and the relative shot-noise (which grows at high fc). 
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Fig. 17. Accuracy of our model and our numerical simulations at redshifts z = 0.35, 1, and 3 for the bispectrum on 
equilateral configurations. We plot our "direct" model (red solid line) and its geometrical modification "tang." (green 
dot-dashed line), shown in Fig. [131 The black dashed line shows the statistical and shot- noise errors of the simulations. 
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Fig. 18. Accuracy of our model and our numerical simulations at redshifts z — 0.35, 1, and 3 for the reduced bispectrum 
on equilateral configurations. We plot our "direct" model (red solid line) and its geometrical modifications "tang." (green 
dot-dashed line) and "flat" (blue dashed line), shown in Fig. [131 The black dashed line shows the statistical and shot- noise 
errors of the simulations. The vertical arrow in the upper right part shows the wavenumber beyond which the simulation 
shot noise is greater than 10%. 
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Combining the halo model used in 
IValaeeas h Nishimichil (j201lt) for the matter power 
spectrum with the 1-loop standard perturbation theory, we 
obtain a good agreement with numerical simulations for the 
bispectrum without any new free parameter. We consider 
the bispectrum as well as the reduced bispectrum, using 
for the latter the power spectrum predicted by the same 
approach (but with the more accurate steepest-descent 
resummation instead of standard perturbation theory). We 
checked that this reasonably good match to simulations 
holds for equilateral as well as isosceles triangles, from 
large to small scales. However, as for the power spectrum, 
the intermediate mildly nonlinear scales are not as well 
reproduced by this direct implementation of our approach. 

The comparison between the three perturbative schemes 
investigated in this article shows that the standard 1-loop 
perturbation theory is actually the most accurate one. 
Because it is also simpler and faster to compute, this is 
also the most efficient one. This is quite different from the 
power spectrum, where the standard 1-loop perturbation 
theory is not the most accurate on large scales and behaves 
badly on small scales (it grows too fast at high fc), so that 
it cannot be used in unified models unless one adds at least 
an external high-fc cutoff. This problem does not occur for 
the bispectrum because the standard 1-loop contribution 
is now negligible on small scales compared with the 1-halo 
contribution. However, if one pushes the perturbative con- 
tribution to higher orders, it may start being too large at 
high k so that one would need to resort to alternative, better 
behaved, perturbative approaches, or to add high-A; cutoffs. 

Next, we have shown how to improve our predictions 
on the transition scales for the matter power spectrum and 
bispectrum, using a simple interpolation scheme instead of 
modifications to halo parameters, which would involve new 
parameters and do not allow significant improvements. This 
method automatically detects this transition regime from 
the shape of the reduced equilateral bispectrum, and in this 
fashion adapts to the change of initial conditions. Thus, for 
CDM linear power spectra that are not pure power laws, the 
transition interval shifts somewhat with redshift 

(i.e., with the local slope n of the linear power spectrum on 
the transition scale) with respect to the interval that would 
be defined by constant thresholds, such as A^(fc±) = 1 
and 30. Moreover, the interpolation through tangent lines 
adapts to the characteristic bends of the power spectrum 
and bispectrum, seen for CDM initial conditions around 
A^(fc) ^ 30. Since this only modifies our model on the 
transition interval [fc_, A;_|_], large scales are still determined 
by systematic perturbation theory and small scales by the 
halo model. Then, we obtain an accuracy on the order of 
10% for the power spectrum and the bispectrum on nonlin- 
ear scales, and 1% on larger weakly nonlinear scales. The 
same levels of improvement and final accuracy are obtained 
for the real-space two-point correlation function. 

Our model can still be improved in various man- 
ners. First, one may investigate other perturbative ap- 
proaches because other resummation schemes may prove 
more accurate than the standard 1-loop perturbation the- 
ory. However, to be more efficient, they should not be 
much more difficult to compute than the standard pertur- 
bation theory. Alternatively, one may go to higher orders. 
For the power spectrum, higher orders are indeed relevant 
because various resummatio n schemes have already been 
shown to be more accurate (jCrocce fc Scoccimarrol [2008t 



iTaruva et ani2009t IValageas fc Nishimichil [201 iD and non- 
perturbative contributions only domin ate after many per- 
turbative orders have become involved (|Valageasll201lh . For 
the bispectrum, the failure of the two resummation schemes 
investigated in this paper to improve over standard pertur- 
bation theory suggests that it may be more difficult to reach 
significant improvements. 

Second, one may improve the halo model used in 
our unified approach. For instance, one could take into 
account substructures (jSheth fc JainI l2003t Giocoli et al.l 
201 0i), devia,tions f rom spherical profiles ( Jing fc Sutoll2002 : 
Sm ith et al.ll2006l) . or the effect of baryons ( Guillet et al.l 
201^^ 

Our model could be extended to other initial conditions, 
especially non-Gaussian one s for which the bispe ctrum is a 
very useful and direct probe (jSefusatti et al.l[2010l) . It would 
also be interesting to use this approach to describe velocity 
fields, and to take into account redshift-space distortions 
([Smith et al.ll2008[) . However, we leave these tasks to future 
studies. 
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Appendix A: Behavior of the bispectrum for one 
low wavenumber in perturbation theory 



5(ki) 




Fig. A.l. A perturbative contribution to the three-point 
connected correlation (5(ki)5(k2)^(k3))c, as in Ea. (IA.4|) . 
The big red circles corresponds to the three nonlinear den- 
sity contrasts 5(ki), (5(k2), and ^(ka). The small black 
dots represent the linear fields (5L(k^) and the Gaussian 
average amounts to connect them by the linear correla- 
tion (5i(k')JL(k;')) = Soi'k'j + \^;)PL{k'j), shown by the 
blue solid lines. There are nn, 71.22, and 71.33 internal lines 
within the circles 5(ki), ^(k2), and (5(k3). There are ni2, 
rii3, and 77.23 lines connecting the three circles. Here we have 
{nil, 7122, ri33} = {2,2,3} and {7112, ni3, 7123} = {3,2,4}. 
This is part of (J^^) (ki)^'") (k2)J(i2) (^3))^, which involves 
the kernels Fg, Fn, and F12. 
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We show in this appendix how the large-scale behavior 
([5^ also arises at all orders of perturbation theory. As is 
well known, within the standard perturbation theory the 
density contrast field (5(k) is written as a perturbative ex- 
pansion over powers of the linear density contrast Sl (k) , 



^(k) = £5(")(k), 



(A.l) 



with 



^(")(k) - y dki..dk„(5,3(ki + .. + k„-k)F„(ki,..,k„) 

x^L(ki).ii(k„). (A.2) 

The kernels Fn can be obtained from recursion rela- 
tions, which are derived from the equatio ns of motion 
(|Goroff et alj 119861: iBernardeau et all l2002al ). and can be 
chosen as symmetric over the wavenumbers {ki,..,k„}. 
Then, the bispectrum can be obtained within perturbation 
theory from the three-point correlation 

(^(ki)^(k2)^(k3))c=5](<5("^nki)^("^Hk2)5("^)(k3))c,(A.3) 



where the subscript c recalls that we only take the con- 
nected part of the Gaussian average. Using Wick's theorem 
and defining the linear power spectrum as in Eq.(IS]) gives 

(^(ki)^(k.)^(k3)). = /ndk5"*^'5^(kf^^+k5"^^^-ki) 



,("13) . ("23) 



(nil) _|,(nil) J^(ni2) ]^("^!') ^ 



("33) _l,(n33) _]^("13) _1j.("23)) 

(A.4) 



where we used a short-hand notation for the diagram shown 
in Fig. lA.ll and the Dirac factors contain sums over the 







x<5i,(-k("" 


)^j^(n23)_ 


X -F2nii+ni2 + 


^l^(«ii) 

"13 V"-j 


X ^2n22+«l2 + 


/1^("22) 
"23 V"-j 


X -P2n33+ni3 + 


,'1.("33) 

"23 y'^j 







wavenumbers k 



(".) 



Of course, the three Dirac factors can 



be combined to factorize out a Dirac prefactor SdO^i + 
k2 -f ka), in agreement with Eq.®. Taking the connected 
part means that each of the three circles in Fig. lA.ll is 
connected to at least one other circle (for instance the case 
nu = Tii3 = n22 = is excluded). 

We are interested in the limit fci — > 0, at fixed ^2 and ^3. 
A well-known property of the kernels f^, which arises from 
momentum conservation (jPeebleslll974t iGoroff et al.lll986l) , 
is that F„(ki -I-.. -|-k„) oc fc^ as k = ki -I-.. -|-k„ goes to zero 
while the individual kj do not, for all n > 2. Then, from 
the first Dirac factor in Ea. (jA.4[) we can see that the first 
kernel F behaves as kf if 2nii -|- ni2 + ^13 > 2. Therefore, 
in the limit fci — ?■ the dominant contributions come from 
the diagrams with 2nii -I- ni2 -I- ni^ — 1, since Fi = 1 and 
for CDM cosmologies Pi,(fci) oc with ~ 1 at low 
ki. This corresponds to either {nu — 0,rii2 = l,"-i3 = 0} 
or {rill = 0, ni2 = 0,ni3 = 1} because we only take the 
connected part of (|A.4p . Thus, for /ci — t- the dominant 



contribution to the perturbative bispectrum reads as 

Spcrt(fcl,fc2,fc3) PUkl) ^I^IYI /n^^"*^ 



X 



("23) 



k2) -P2"33-|-n23 



^J^("33) _J^("33) _Jj;("23)-j 



p ,'1^("22) . ("22) 1 u("23)n 

^2"22 + l + "23lJ^J ^-'^J '^^I'S' 



n^i(fci"*') + lpcrm. 



where we note by "1 perm" the 
tion with respect to {k2 ks}. 
in Eg. ljA.Sp . because the limit at 

("22) i,(n22) 



(A.5) 



symmetric contribu- 
There is a subtlety 
low fci of the ker- 

,("23)^ 



nel F2„,,+i+„,3 -k;"^^^ -ki , k) "^") contains diver- 

gent terms of the forn i (ki ■ k^ ) / k^ (jGoroff et al.l 119861: 
IScoccimarro fc FriemanI [l996[ ). However, taking the limit 
fci — >■ elsewhere, we obtain k3 = — k2 so that the symmet- 
ric term with respect to {k2 O ks} actually corresponds to 
a change of sign of all wavenumbers. This cancels the di- 
vergences of the form (ki • ]s.j)/kf so that the limit fci — 
is finite. The fact that such infrared divergences cancel out 
for equal-time statistics can be trace d to the Galilean in- 
variance of the equa tions of motion (jJain fc Bertschingerl 
119961: IValagea^l2n04ll . 

Terms at all orders of perturbation theory contribute 
to the partial large-scale limit (jA.5|) through the nonlin- 
ear corrections to ^(k2) and ^(ka). In contrast, as we have 
seen above, only the linear term SlO^i) contributes. The 
behavior 



(IA.5I) can also be understood from the physical 
arguments discussed in Sect. 12.31 and the explicit result 
(|A.5P confirms that this asymptotic behavior should be 
preserved by both high-order perturbative terms and non- 
perturbative corrections. 

It is interesting to note that the "renormalization" of 
the prefactor to the -P(fci) tail by the higher order terms 
in Ea. (jA.5|) never occurs for the power spectrum P{k), 
where we recover P{k) Phi k) at low k. Th ere, as seen 
from the analysis described in IValaeeaj (|2002f ) , perturba- 
tive terms beyond linear order scale at least as k^Phik) 
at low k. This can also be understood from the fact that 
for the power spectrum there are no partial large-scale lim- 
its. On the other hand, if all wavenumbers go to zero, as 
fci ~ fc2 ~ ^3 ~ with fc — > 0, higher order contributions to 
the bispectrum scale at least as fc^PL(fc)^ while the "tree- 
order" result scales as Pi(fc)^, so that in the simultaneous 
large-scale limit we also recover the lowest order result. 

Appendix B: Two-point correlation function 

Although the main focus of this paper is the computation 
of the bispectrum, we have seen in Sects. [7] and [S] that the 
shape of the reduced bispectrum can be used to improve the 
model devised for the power spectrum. This leads in turn to 
an improved model for the real-space two-point correlation 
function, which is given by 

sin(fca;) 



f (x) = 47r / dfc fc^ P{k) ■ 



fca 



dfc , 9 , sm{kx) 
k kx 



(B.l) 
(B.2) 



Therefore, we compare in this appendix the two-point cor- 
relation functions obtained either with our "direct" model. 
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Fig. B.3. Accuracy of our models and our numerical simulations at redshifts z ~ 0.35, 1, and 3 for the real-space two- 
point correlation function. As in Figs. IB.l] and [R2l the red solid line "direct" is the relative difference P5|) between the 
simulations and our model as described in , Valageas fc Nishimichi (2011|1. while the green dashed line "tang." corresponds 
to the geometrical modification shown in Fig. [M) The black dashed line shows the relative statistical error (which grows 
at large x) and shot-noise error (which grows at small x) of the simulations. 




0.1 1 10 60 90 100 110 120 

X [hi-i Mpc] X [hi-i Mpc] 



Fig. B.l. Real-space two-point correlation function ^(x), 
at redshifts z = 0.35, 1, and 3 (from top to bottom). We 
show the "direct" model (red solid line) and its geometri- 
cal modification "tang" (green dot-dashed line). They are 
the Fourier transforms of the corresponding power spectra 
shown in Fig. [TH 



Fig. B.2. Real-space two-point correlation function ^(x) 
at redshifts z = 0.35, 1, and 3 (from top to bottom) as 
in Fig. IB. II but on larger scales. In addition to the "di- 
rect" model (red solid line) and its geometrical modifica- 
tion "tang" (green dot-dashed line), we also show the linear 
two-point correlation (black solid line). 



already studied in 'Valag eas &: Nishimichil (|201lf ) , or with 
the geometrical modification "tang" to the power spectrum 
explained in Fig. [HI Thus, we plot in Figs. lB.ll and lB.2l the 
two-point correlation functions defined by the two power 
spectra of Fig. [HI through the Fourier transform (|B.2p . As 
expected, we can see in Fig. IB. II that, as for the power 
A^(A;) in Fig. [Ml the geometrical modification "tang" cor- 
rects most of the artificial "dip" that was shown by the 
"direct" model on transition scales (i.e., the mildly non- 
linear regime). The curves even look closer to the simula- 
tion results and more regular because the integral (IB.2P has 
regularized the geometrical modification. Indeed, while the 
power Aj3^j^g(fc) shown in Fig. I14l was only continuous, with 
a discontinuous derivative at the correlation ^tang(2;) 
has a well-defined and finite derivative over all cc > 0. 



Contrary to the power spectrum, the modification 
"tang" to the two-point correlation is not restricted to a 
finite range [x-, x+\ because the modification of the power 
on the wavenumber interval [fc_ , fc+] contributes to all dis- 
tances X through the integral (|B.2I) . However, as expected, 
we can check in Fig. IB. II that on very small and very large 
scales the modified correlation f tang (2;) becomes increas- 
ingly close to the original model prediction ^direct (a^), and 
as such shows a good agreement with numerical simula- 
tions. 

Nevertheless, to check in more detail that the modifi- 
cation "tang" does not spread on large scales, where the 
initial power is much weaker, we show in Fig. IB. 21 the cor- 
relation functions of Fig. IB. II on larger scales, around the 
baryon acoustic oscillation. We can see that both corre- 
lations, '^tang(a^) and ^direct (a^), are very close and within 
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the error bars of the numerical simulations. In particular, 
they capture very well the departure from the linear corre- 
lation, mostly due to the 1-loop perturbative contribution. 
Therefore, the modified correlation ^tang(a;) provides a good 
model from small to very large scales. 

As in Fig. [12] for the power spectrum, we show in 
Fig. IB. 31 the relative accuracy of our analytical models and 
numerical simulations. (The curve obtained for the "di- 
rect" model was already shown in 'Valage as fc Nishimichil 
(|201lt) .) In agreement with Figs. [TBI and lB.li we can see that 
the modification "tang" provides a significant improvement 
on transition scales. Thus, it ensures an accuracy of 10% or 
better on nonlinear scales, while an accuracy on the order 
of 1% is again reached on quasilinear scales, associated for 
instance with the baryon acoustic oscillation of Fig. IB. 21 
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